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ABSTRACT 


A  numerical  integration  procedure  which  employs  m  values  (ordinates) 
of  the  integrand  is  known  as  an  m-point  formula.  More  generally  if  the 
integrand  is  w(x  )  f(x  ,...,x  )  where  w( • )  is  called  the  weight 

function,  we  could  have  an  m-point  formula  employing  m  values  of  f(*) 
to  approximate  the  integral  over  a  region  in  n-space.  The  m-point  formula 
is  said  to  have  precision  p  if  it  gives  exact  results  whenever  f(*)  is 
a  polynomial  of  degree  <  p.  One  dimensional  integration  formulas  can  be 
used  to  obtain  a  formula  over  a  rectangular  region  in  n-space  by  forming 
products  of  the  integrals  over  each  separate  variable. 

For  example,  a  Gaussian  m-point  formula  in  one  variable  achieves  a 
precision  2m-l.  We  can  obtain  from  this  an  m’^-point  formula  of  precision 
2m-l  in  n  dimensions.  Although  this  formula  is  not  unduly  extravagant 
in  points  for  fixed  n  and  large  p,  one  soon  has  an  impractical  problem 
if  the  number  of  dimensions  n  becomes  large.  A  formula  economical  in 
points  for  fixed  p  and  large  n  is  thus  of  considerable  practical  value. 

The  more  important  known  results  in  this  field  are  as  follows: 

1)  An  n  +  1  -  point  formula  of  precision  2  valid  over  an  arbitrary 
region  in  n-space--weight  function:  arbitrary. 

2)  An  n  +  2  -  point  formula  of  precision  5  the  n- simp lex- -weight 
function:  1. 

5)  A  2n-point  formula  of  precision  3  for  an  arbitrary  symmetric  region 
is  n-space- -weight  function:  symmetric  but  otherwise  arbitrary. 

4)  A  2n^  +  1  -  point  formula  of  precision  5  the  n-cube  and  the 
n-sphere--welght  function:  1. 

5)  Results  on  transformations  and  symmetric  regions. 
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The  following  results  presented  in  this  thesis  are  believed  to  be  new. 

1)  A  2n^  +1  -  point  formula  of  precision  5  f<5r  an  arbitrary  fully 

symmetric  region  in  n-spacc. 

2)  Two  "Ipn^  +I4n^  "T^i)  +1  "  point  formulas  of  precision  9  valid 

over  an  arbitrary  fully  symmetric  region  in  n-space.  For  each  of 
these  formulas  the  weight  function  is  fully  symmetric  but  otherwise 
arbitrary.  Results  of  four  particular  regions  and  weight  functions 
are  tabulated  to  25  dimensions  in  the  Appendix. 

These  integration  formulas  were  obtained  by  solving  a  system  of  non- 
linear  algebraic  equations.  Polynomials  in  one  variable  orthogonal  over  the 
fully  symmetric  region  with  respect  to  the  fully  symmetric  weight  function 
were  an  aid  in  obtaining  the  integration  formulas  since  the  zeros  of  these 
polynomials  can  be  shown  to  be  the  non-linear  unknowns  in  the  non-linear 
algebraic  equation. 

For  large  n,  the  method  of  attack  developed  here  will  in  general  lead 

2k 

to  integration  formulas  of  precision  p=  l<^k  +  1  (k  an  integer)  using  0(n  ) 

points j  that  is,  the  formulas  are  economical  for  large  n  and  fixed  p. 

3)  A  -point  formula  of  arbitrary  precision  is  also  developed 

here  for  the  finite  and  infinite  n-sphere. 

4)  Estimates  of  the  error  committed  in  performing  numerical  inte¬ 
gration  are  eminently  desirable  but  not  always  easily  obtained. 
Towards  this  end  contour  integration  and  asymptotic  techniques  were 
employed  to  extend  known  results  of  error  bounds  for  one-dimen¬ 
sional  integration  to  error  bounds  for  repeated  Gaussian  integration 
in  the  case  when  the  integrand  is  a  meromorphic  function  of  n 
complex  variables. 

A  number  of  examples  are  given  throughout  the  thesis  and  in  the 
Appendix  illustrating  the  accuracy  of  the  formulas  developed. 
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CHAPTER  1 
INTRODUCTION 


This  introductory  chapter  is  intended  to  facilitate  the  understand¬ 
ing  of  the  thesis. 

Most  integrals  cannot  be  evaluated  in  terms  of  simple  functions. 

This  difficulty  is  discussed  in  section  l)  and  points  to  the  need  for  numerical 
integration  formulas. 

The  second  section  contains  (i)  some  basic  definitions  of  terms  used 
throughout  the  thesis,  and  (ii)  some  examples  of  numerical  integration  formulas. 

Although  the  study  of  numerical  integration  formulas  has  received 
considerable  attention,  many  difficult  problems  still  remain  unsolved^  this 
being  particularly  so  for  higher  dimensions.  Section  3)  lists  the  problems 
this  thesis  is  concerned  with. 


1)  The  Purpose  of  Numerical  Integration  Formulas 


The  evaluation  of  integrals  often  poses  a  difficult  problem.  Although 
we  can  evaluate  a  lot  of  integrals  in  terms  of  simple  functions j,  most  of  the 
integrals  arising  in  research  or  industry  cannot  be  evaluated  in  terms  of 
simple  functions,  and  can  hope  to  evaluate  them  only  using  a  numerical 
integration  formula.  This  is  especially  true  of  multidimensional  integrals, 
due  to  the  multidimensional  complexity  of  the  integrand. 


For  example,  in  one  dimension  we  have 
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and  80  on,  all  of  these  being  Integrals  that  cannot  be  evaluated  in  terms  of 
simple  functions.  Note  that  we  can  differentiate  every  one  of  the  intiO'giands 


2  - 


as  often  as  we  please,  since  continued  differentiation  does  not  give  rise  to 
new  transcendents.  Continued  integration,  however,  soon  does.  Hence  if  a 
new  integral  arises  in,  say,  research,  one  is  lucky  to  be  able  to  evaluate 
it  in  terms  of  simple  functions  in  one  dimension,  and  if  it  is  multidimensional, 
he  is  so  much  more  likely  to  run  into  difficulty. 

2)  Definitions  and  Examples 

i)  Definitions.  In  order  to  facilitate  further  discussion,  we  now 
make  the  following  definitions. 


(2) 


dX 


n 

That  is,  X  is  a  1  x  n  row  vector  in  euclidean  n-space  E  ,  dX  is  a  hyper= 


volume  element  in  this  space. 


We  want  to  find  numerical  integration  formulas  of  the  form 


(3) 


R 

n 


m 


Here  f(x^ . x*^)  Is  the  continuous  function  we  wish  to  Integrate  over  a 

region  In  with  respect  to  the  continuous  weight  function 

w(x^,...,x^)  that  does  not  change  sign  over  the  region. 
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In  terms  of  (2)  we  may  write  (j)  In  the  foim 
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Formula  (4)  is  said  to  be  of  precision  p  if  it  gives  exact 
results  whenever  f(X)  is  a  polynomial  of  degree  <  p.*  That  is,  the  gener¬ 
al  term  of  such  a  polynomial  may  be  written 


(5) 


n  -t  ^ '  Y"’ 

H  (x  )  ^  ,  where  0  <  y  ^  P  >  - 

i=l 


0 


When  we  are  looking  for  integration  formulas  of  the  type  (4),  we 

have  the  problem  of  finding  the  constants  c^  which  should  preferably  be 

real,  and  the  points  which  should  also  be  real  and  lie  within  the  region 

R  .  Moreover,  given  p,  we  want  to  minimize  m  . 
n 

ii)  Examples .  Some  examples  of  numerical  integration  formulas  are 
Gaus  s -Legendre „ 


Here  w(x)  =  1,  p  =  2m“l,  c^  =  2/{  (l-x^^)  [P^(x^ )  ]^) ,  the  ^ j  '  ® 
m  zeros  of  the  polynomial  orthogonal  over  the  interval  (-1,1) 

with  respect  to  the  weight  function  1,  As  an  illustration  we  have 


(T)  J  f(x)dx  ^  I  +  9  +  9 

-1 


Chebychev  (one  dimension).  The  following  equi-weighted  formula  is  due  to 
Chebychev  [25]. 


-1  j  =  l 


* 


We  will  later  show  that  this  definition  does  not  "fit"  very  well. 
A  discussion  of  precision  is  given  in  Chapter  III. 
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Here  w(x)  =  l/sTl  -  ,  p  =  2m  -  1  ,  *^he  ™  zeros 

of  T^(x)  orthogonal  over  the  interval  (-1,1)  with  respect  to  the  weight 
function  l/Jl  -  x^  .  As  a  more  specific  example  we  have 


(9) 


f  ^  5  (£(-|/f)  +  £(o)  +  £(7t)l 

-1 


where  m=3,  p=5,  c=c=c 


1  -  -2  -  -5  =  3-  \  =  -Jk'  ^2  =  =^5  =  '/F  • 


Tyler’s  formula  for  the  n-cube. 

(10) 


,n-l 


2n 


r  '‘-I  f(x^, . . .  ,x’^)dx^  ...  dx’^  =  ■—“  \  f(x^,o.,,x”) 

J  .1  J  .1  J 

j  =  l 


Here  w(X)  =  1,  m  =  2n,  p  =  3,  c  =  x^  =  /| 
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6,  .  being*  we  Unknown  Kronecker  delta. 
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3)  The  Problem. 


One  way  of  obtaining  a  numerical  integration  formula  over  a  rectang^ 
ular  region  in  n  dimensions  is  by  taking  products  of  one  diiiiensional  integ= 
rals  over  each  separate  variable.  Proceeding  thus,  we  would  find  the  number 
of  integration  points  increasing  exponentially  in  n.  For  example,  by  taking 
products  of  formulas  of  the  type  (6)  and/or  (8)  and  keeping  the  precision 

Hi 

p  =  2m-l  in  each  variable,  we  would  require  m  points  over  a  rectangular 
region  in  n  dimensions.  As  the  number  of  dimensions  Increases  we  would 
soon  find  ourselves  faced  with  a  problem  that  takes  too  long  to  evaluate 
even  with  the  fastest  digital  computer.  Hence  there  is  a  need  for  numerical 
integration  formulas  that  require  fewer  points  for  a  large  number  of  dimens^ 
ions  n. 


Numerical  integration  formulas  of  precision  p  in  n  dimensions 
requiring  fewer  than  (£±1)"  points  have  been  found  by  earlier  workers  in 
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the  field  but  only  up  to  precision  three  in  the  general  case  and  up  to 
precision  five  in  two  particular  cases.  Such  a  limited  class  of  formulas 
severely  restricts  the  accuracy  obtainable  when  performing  numerical 
integration. 


Although 


point  formulas  require  an  extremely  large 


number  of  points  when  the  dimension  n  is  large,  we  cannot  overlook  the 
superior  accuracy  of  these  formulas  for  a  large  class  of  functions  that 
have  all  their  derivatives  continuous.  Moreover,  these  are  presently  the 
only  available  formulas  of  arbitrary  precision.  Of  the  many  different 
regions  possible  in  higher  dimensions,  we  have  seen  only  such  formulas 
(of  arbitrary  precision)  for  rectangular  regions  in  n-space. 


Very  little  work  has  been  done  regarding  error  estimates  of 


integration  formulas  in  n  dimensions.  This  may  be  expected  since  the 
topic  is  relatively  new.  Most  workers  have  concentrated  on  the  task  of 
finding  workable  results  and  in  general  the  results  known  are  of  such 
complexity  that  practical  error  evaluation  would  be  difficult  to  obtain. 
Apriori  error  estimates  are  preferable  when  they  can  be  found,  particularly 
if  the  time  required  to  evaluate  these  is  small  compared  with  the  time 
required  to  evaluate  the  approximate  integral. 
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CHAPTER  II 

REVIEW  OF  THE  LITERATURE 


This  chapter  of  the  thesis  states  the  known  results  of  the.  field,, 
It  is  by  no  means  a  complete  summary  of  approximate  integration,  but  c.on“ 
sists  rather  of  known  results  of  n-dimensional  approximate  integration  per^- 
taining  to  the  problems  discussed  in  Chapter  I,  Stroud  [24]  gives  a  very 
extensive  bibliography  of  approximate  integration  up  to  I96O. 

In  section  l)  we  state  results  of  nineteenth  century  authors, 

C,  F.  Gauss,  P.  L.  Chebychev  and  J.  C.  Maxwell, 

The  principal  results  of  section  2)  are  concerned  with  (a)  ex“ 
tending  integration  formulas  known  over  one  region  to  other  regions,  (b) 
construction  of  integration  formulas  for  basic  regions  such  as  n-cubes, 
n-simplexes,  n-spheres,  (c)  the  construction  of  integration  formulas  for 
arbitrary  regions  in  n-space,  (d)  numerical  integration  by  Monte  Carlo 
Method,  and  (e)  a  derivative  error  bound  of  von  Mises. 


1)  Nineteenth  Century  Results 

Gauss  [1]  was  the  first  to  solve  the  problem  of  obtaining  maximum 
possible  accuracy  with  m  points  in  one  dimension.  By  chosing  his  m 
points  Xj  to  be  the  m  zeros  of  P^(x),  he  was  able  to  obtain  precision 
p  =  2m  -  1. 

J.  C.  Maxwell  [2],  in  I877,  found  numerical  integration  formulas 
for  the  rectangle  and  for  the  5-dimensional  parallelelepipedon.  For  the 
rectangle,  he  obtained  a  formula  of  precision  7  using  I5  points; 


(1)  J"  f(x,y)dx  dx  =  P  f(0,0)  +  f(+p,0)  +  f(0,+p)  +  R  f(±qs±r) 


-1  -1 


(-1.1)  (1,1) 


In  order  to  set  up  equations  and  solve 
for  the  points  p,  q,  r  and  the  numbers 
P,  Q,  R  avC  shown  in  the  figure  on  the 
left,  he  first  transformed  the  rectangle 
into  the  square. 

He  then  found  the  solution  to  the 
following  six  simultaneous  non-linear 
algebraic  equations. 


1 


i 


} 
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P  +  4Q  +  8R  =  U 

Qp^+2R(q2+r2)  =  2/3 

^2)  Qp^2R(qVr^)  =  2/5 

Qp^+2R(q^+r^)  =  2/7 

2R  q^  =  I/9 

Rq2r2(q2+r2)=  1/15 


The  solutions  he  obtained  are 


p2  =  12/55 

(5)  =  3/5  [1+(6/5iF] 

=  3/5  [1“( 6/31)2] 


p  =  8/162 
Q  =  98/162 

R  =  31/162 


For  the  three  dimensional  case  of  precision  7  he  obtained  two  27 
point  solutions,  each  of  which  had  points  that  lay  outside  the  region  of 
integrationo  This  sort  of  thing  occurs  quite  frequently  when  one  tries 
to  find  integration  formulas  in  more  than  two  or  three  dimensions  (see  for 
example,  equation  (lO)  of  Chapter  I  for  n  >  5)»  By  a  different  selection 
of  points.  Maxwell  could  have  avoided  obtaining  points  that  lie  outside 
the  region  of  integration  (see  M.T.A.C.  v.  12  P.  27I+).  Whenever  possible 
we  want  to  avoid  obtaining  formulas  with  such  points  since,  for  example, 
if  we  wanted  to  compute  the  weight  of  a  petroleum  product  by  measuring 
its  density  which  varied  throughout  a  rectangular  tank,  the  formulas  command 
us  to  measure  the  density  outside  of  the  tank. 


For  the  same  reason,  we  want  to  try  to  avoid  obtaining  formulas 
with  complex  points.  An  example  of  this  phenomenon  is  furnished  by  the 
Chebychev  type  formula; 


(An  integration  formula  such  as  this,  for  which  all  the  c.’s  are  equal, 

J 

is  said  to  be  of  a  Chebychev  type,  in  honor  of  Chebychev  who  first  tried 
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to  find  such  integration  formulas).  For  the  integration  formula  (l4)  to 
have  precision  of  at  least  m  for  the  particular  case  w(x)  =  1,  the 
Xj ' s  can  be  shovm  (see  [3])  to  be  the  m  zeros  of  the  polynomial  part  of 


(5) 


exp 


ra 

2 


J 


r  iog( 


X 


t  )dt 


We  will  call  the  polynomial  part  of  (5)  G^(x)„  Chebychev  found  that  all 

the  polynomials  G^,  G^,  G^  and  G^  had  real  zeros  within  the  inter- 

val  (-1,1).  Gg,  however,  had  six  complex  zeros.  Chebychev  stopped  here, 
leaving  the  question  as  to  what  the  zeros  of  ^t’e  like  for  m  >  10 

open.  Later,  after  Chebychev' s  death  it  was  shown  [4]  that  for  m  >  10 
every  ^t  least  one  pair  of  complex  roots. 


2)  Recent  Developments 

The  problem  of  obtaining  nximerical  integration  formulas  in  higher 
dimensions  received  little  attention  until  the  last  20  years.  Indeed  no 
practical  use  could  be  found  for  these  until  the  development  of  electronic 
computers  (around  1947)  made  it  possible  to  envisage  extensive  computations 
of  integrals  in  more  than  one  variable.  Since  then  new  attempts  at  finding 
numerical  integration  formulas  in  higher  dimensions  have  been  made. 


The  following  definitions  facilitate  the  presentation  of  the 
results  of  this  section. 

Definitions 

(a)  A  set  R  in  is  said  to  be  fully  symmetric  if  X  e  R  implies 

Y  e  R  where  Y  is  any  point  obtainable  from  X  by  permutations  and  by 
n 

changes  of  sign  of  the  coordinates  of  X. 


(b)  A  function  g  defined  in  a  fully  symmetric  set  is  fully  symmetric  if 

g(x)  =  g(K). 
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(c)  A  numerical  integration  formula  is  said  to  be.  fully  symmetric  if  the 

sets  of  evaluation  points  X.,  X  ,  X  tj.v'e  decomposable  into  fully 

12  m 

symmetric  sets  and  if  the  weight  function  w(X)  is  a  fully  symmetric  function. 

Tyler  [5]  in  1953  found  integration  formulas  over  regions  bounded 
by  parabolas  and  the  hypercube.  His  numerical  integration  formula  for  the 
hypercube  is  of  precision  5,  very  similar  to  formula  (lO)  of  Chapter  I,  and 
requires  2n  +  1  points. 

Hammer,  Marlowe,  Stroud  and  Wymore  (but  mainly  Hammer  and  Stroud) 
have  written  extensively  on  the  problem  of  finding  numerical  integration 
formulas  in  higher  dimensions.  They  have  given  numerical  integration  form¬ 
ulas  up  to  precision  3  for  the  n-simplex  [6]  9  up  to  precision  5  foi^  the  n- 
cube  and  the  n-sphere  [7];  they  have  given  several  theorems  for  extending 
available  results  into  higher  dimensions  and  into  regions  related  to  given 
regions  by  transformations.  We  present  here  two  important  theorems  due 
to  Hammer  and  Wymore  [8]. 

i)  Transformation  theorem 

Suppose  we  have  an  integration  formula  which  may  be  written  in 
terms  of  its  error 

m 

(6)  E(£.R„)  =  y  -  f  "it*)  ® 

j=l  '  \ 

where  w,  and  f  are  continuous  functions  of  X  in  a  region  of 

1  n 

euclidean  n-space  e”. 

Let  S  be  another  region  in  E*'  such  that  there  exists  a  trans- 
n 

formation 

Y  .  Y(X)  (.  y^(X),  y"(X)) 


(7) 


j,  iaoxtyj'Ms.frt  A  (o) 


n/>,i 


.  Ij 


''-‘P.H 


>!•  ^.17J  -LrtW'ryi 


‘T 


M-rv:: 


,i 


:'  ij'.' 


r 


0  \:U’  : 


•'  ’  -C  -> 


':.  f 


0 


"■  1  O  ^.;j 

,,’V  - 

'  & 

C<  ■ . 

# 


«Oi 


10 


with  a  continuous  non-vanishing  Jacobian 


(8) 


J  =  J(x)  = 


hi 

3x1 


Sx" 


Ox 


which  maps  R  onto  S 
n  n 


Further,  if  w^(X)  =  W2(Y)  and  g(Y)  is  a  continuous  function  of 


Y  in  S  ,  we  have 
n 


(9) 


f  w,(^)  g(V)  dY  =  r  w  (Y)  g(Y)  |J(X)|  dX 
dS  d^ 


n 


n 


=  f  w  (X)  h(X)  lJ(X)|  dX 

Jr  ^ 


n 


where  we  have  let  g(Y)  =  h(X). 

Hence,  by  formula  (6)  we  have 
m 

(10)  E(|Jlh,R^)  =  ^  aj|j(Xj)|  -  Wg(Y)  g(Y)  dY  . 

j=l  ‘ 

By  formulas  (lO)  and  (6)  we  thus  have  the  following 

theorem; 

m 

(11)  E(g.S„)  =  /  b,  g(V  )  -  f  Wg(Y)  g(Y) 


i  S'  Js 

j«l  n 


dY 


where  bj  a  aj|j(Xj)l,  Yj  a  hence  for  every  formula  of  the  type 

(6)  in  R  there  is  a  corresponding  formula  of  the  same  type  in  S  with 


error  function  E(g,S^)  «  E(|j|h,R^) 


i 


11 


COROLLARY:  If  Y(X)  is  a  non-singular  transformation  for  which  j(X)  is  a 

constant,  then  E(g,S^)  =  E(|jjh,R^)  =  jj]  E(h,R^)  and  hence  if  formula  (6) 

has  precision  p  in  R  ,  then  formula  (ll)  has  precision  p  in  S  .  In 

n  n 

particular,  if  Y(X)  is  an  affine  transformation,  J(X)  is  a  constant. 


ii)  Cartesian  product  theorem 

Let  R  be  the  Cartesian  product  of  two  regions  R^  and  R^  in 
lower  dimensional  euclidean  spaces.  Let  us  designate  X  in  R^,  Y  in 
R  ,  so  that  every  Z  in  R  can  be  written  (X,Y).  Suppose  we  have  numer” 
leal  integration  formulas  in  R^  and  R^  with  weight  functions  w^(X)  and 
w^Cy)  respectively.  Together  with  their  error  functions,  these  may  be 


written 


m. 


(12)  °  Z  ® 

1=1  '^’^1 


m. 


(13) 


J  2  j  Jr  2  2 

j=l  2 


Using  (12)  and  (I5)  we  obtain  the  following 
THEOREM:  Let  f(z)  =  f(Xs,Y)  be  defined  and  continuous  on  the  region 

R  =  Rj^  X  Rg,  Then,  if  we  define 

Tti.  m 
1  id 

(11»)  E(R,£)  =  y  ^  a.  bj  £(X..Yj)  -  f  J  w^lx)  Wg(Y)  f(X,Y)  dX  dY 


1=1  1=1 


R^  ^R, 
2  1 


we  have 


m. 


(15)  E(R,£)  -  ^  bj  E^(R^,£(X,Yj))  +  f  w^(X)  Eg(Eg.£(X,Y))  dX 


J=1 
m. 


Z  ”1  f  Ej(8i,£(X,Y))  dY  . 


’)0 


V 


'-'f- 


rW 
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COROLLARY  (l)!  If  f(X,Y)  is  a  function  such  that  E^(R^,  f  (X^,Y))  =0  for 
each  Y  e  R^  and  E2(R2>  f  (X,Yj ))  =0  for  each  X  e  R^^  then  E(R,f)  =  0 
and  formula  (lU)  is  exact  for  f(X,Y). 


Now,  if  is  a  class  of  functions  defined  on  R^  and  is  a 

class  of  functions  defined  on  R^  then  the  Cartesian  product  class 

F  =  F.  X  F  of  functions  defined  on  R,  x  R  is  the  class  of  all  functions 
12  12 

f(X,Y)  such  that  f(X,Y)  e  F^  for  each  Y  in  R^  and  f(X,Y)  e  F^  for 

each  X  €  R  . 

1 


COROLLARY  (2):  If  classes  of  functions  F  and  G  respectively  defined 


on 


R^  and  R^  are  representable  as  all  linear  combinations  of  basis  sets  of 

functions  f^,.,.,f^,  respectively  then  (F  x  G)  is  the  set  of 

all  functions  with  f,  g.  as  basis. 

i  j 

PROOF  (of  corollary  (2)): 

r  s 

If  h  -  c,,  f,  g.  where  the  c,.  s  are  constants,  then 

i=l  j=i  ij  i  ij 

h  e  F  X  G.  To  show  that  h  €  F  x  G  implies  h  is  a  linear  combination  of 

r  s 

the  (f.g.)i.  we  observe  that  h  =  Z  a.f.  =  .Z  b,g.  where  the  a 

1  j  i=i  iij— ijj  1 

unique  functions  on  and  the  b^  are  unique  functions  on  R^. 


are 


Let  X  , .,.,X  be  the  points  in  R-  such  that  the  determinant 
Is  1 

|f^(Xj^)j  ^  0.  Then  we  have  with  b^^  the  value  of  b^  at  X^, 
r  s 

^Z,  a^f^(X,  )  =  .Z,  b.,  g.f  k  =  which  are  identities  on  R^.  From 

ial  1  i  k  jal  JK  J  2 


this  set  of  equations  we  can  solve  for  a^  b  ^Z^^  c^^j  gy  the  c^^ 
unique  sine©  the  determinant  ^  • 


s  being 


Hi)  Implications  of  the  theorems 

What  are  the  implications  of  the  above  theorems?  The  transforma¬ 
tion  theorem,  tells  us  that,  for  example,  if  we  have  a  numerical  integration 
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formula  of  precision  p  over  a  given  region,  we  automatically  have  an  inte¬ 
gration  formula  of  the  same  precision  over  a  region  related  to  the  given  one 
by  an  affine  transformation.  This  result  is  not  only  useful  in  extending 
known  integration  formulas  to  regions  related  to  a  given  one  by  affine  trans¬ 
formations,  but  it  is  also  a  tremendous  aid  in  obtaining  new  integration 
formulas . 


To  illustrate  this,  the  number  of  different  monomials  in  the 


polynomial 


(1  + 


1  2 
X  +  X  + 


n\p 

+  X 


is  equal  to  the  number  of  ways  of  assigning  p  indistinguishable  objects 
to  (n  +  l)  different  groups,  which  is  well  known  to  be  the  binomial 
coefficient , 

+  p^^  (n+l)p  (p+l)^  ^  + 

py  p!  "  n!  "(^n 

Hence  finding  an  integration  formula  of  precision  p  over  an  arbitrary 
region  in  would  require  the  solution  of 

(n+l) 

_ _ E 

pI 

simultaneous  non-linear  algebraic  equations.  This  number  of  equations  not 
only  becomes  very  large  as  n  or  p  increases,  but  it  is  dependent  on  n 
and  p. 

If,  however,  we  restrict  ourselves  to  a  fully  symmetric  region 
when  looking  for  integration  formulas,  relying  on  the  above  transformation 
theorem  to  extend  our  results  to  other  regions,  the  total  mimber  of  simult¬ 
aneous  non-linear  algebraic  equations  becomes  independent  of  n,  and  con¬ 
siderably  less.  For  example,  to  obtain  an  integration  formula  of  precision 


i  ;■'  •  :■  j 


j. 


i  i_ 


r  ')■ 


.  ; 


Q;r:v.: 


'x.-'o  6):]. *4 ?“/:'«• 

'  .Beet 
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7  over  arbitrary  5-space  would  require  the  solution  of  120  simultaneous  non¬ 
linear  algebraic  equations.  For  a  fully  symmetric  region  however,  (in  3- 
space  as  in  n-space)  the  total  number  of  simultaneous  non-linear  algebraic 
equations  reduces  to  7, 


The  second  theorem  tells  us  that  if  we  have  an  integration  formula 
over  the  circle  and  another  over  an  interval  along  a  line,  we  immediately 
have  one  for  the  cylinder;  other  examples  could  easily  be  adduced.  To  obtain 
an  integration  formula  over  the  square,  we  would  simply  use  the  fact  that  the 
square  is  a  Cartesian  product  of  two  lines  of  the  same  lenth,  together  with 
a  Gauss-Legendre  integration  formula  along  a  line. 

Note  (1);  If  we  need  m^  points  to  obtain  precision  p^^  in  the  variables 

X  in  R  and  m  points  to  obtain  precision  p  in  the  variables  Y  in 

12  2 

R  ,  then  we  need  m.  x  m.,  points  to  obtain  precision  p,  in  the  variables 
2  12  i 

X  and  precision  p^  in  the  variables  Y  in  the  region  R^  x  R^.  That  is, 
the  precision  we  can  guarantee  in  our  Cartesian  product  integration  formula 
is  the  minimum  of  p^  and  p^. 

iv)  A  result  for  cones 

We  now  develop  an  interesting  integration  formula  for  cones,  due 
to  Hammer,  Marlowe  and  Stroud  [6]. 


Let  R  be  an  n-dimensional  region  on  the  hyperplane  y  =  1  which 
n 

is  embedded  in  We  represent  the  points  in  by  (X,y),  X 

being  a  point  in  E^.  Then  the  set  of  all  points  yR  ,  where  0  <  y  <  1 

n+l 

is  a  cone  C  ,  with  base  R  and  vertex  at  the  origin  in  E  .  Let 
n+1  n 

f(X,y)  be  a  function  defined  over  and  suppose  we  have  a  numerical 

integration  formula  over  the  base  R  of  C  , 

°  n  nf  i 


If,  for  example 
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m. 


(16) 


then 


(it) 


J  f(X,l)  dX  ^  ^  f(X^,l) 

\  i=l 

f  f(X,y)  dX  dy  =  r  f  f(X,y)  dX  dy 

,  J  0  ^  yR 

n+1  ■'  n 


m. 


'  /  I 


dy 


0  i=l 


since  the  Jacobian  of  the  affine  transformation  from  R  to  yR  is  y 

n  n  ^ 


■>0 


If  we  now  define 


m. 


(18) 


we  have 


(19) 


g(y)  =  ^  f(Xiy»y)  » 

i=l 

f  f(X,y)  dX  dy  =  f  y”  g(y)  dy  . 
n+1 

Since  numerical  integration  formulas  of  the  form 
1 

(20)  I'  y"  g(y)  dy  =  ^  g(yj) 

'  0  j.l 

certainly  exist,  we  obtain  the  result 

”*1  ^2 

(21)  J  £(X.y)  dX  dy  2r  ^  X  “l  ‘’j  ' 


n+1 


isl  Jal 


Hence  if  (16)  holds  precisely  for  polynomial  f(X)  of  n  vari¬ 
ables  of  at  most  degree  p  over  a  region  R  ,  and  if  formula  (20)  holds 
precisely  for  polynomials  g(y)  of  at  most  degree  p  in  y,  then  (21) 
holds  over  the  cone  for  all  polynomials  containing  monomials  of  at 

most  degree  p  in  its  n+1  variables. 


.1, 


„  i 

j  ' 


(rr; 


i, 


i  ■ : 


t  ;•;, 


'  s  '*  'iS'- ^  tr-siwo  .’ij ,  :3i*  VwJI'Ij# 

......lA.,  XIa  tot  J),  sd;j  t®vc  a&lori 
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v)  A  result  for  simplexes 

n  n 

Let  the  vertices  of  the  n-slmplex,  S  be  Y, ,  , . . ,  Y  , „  and 

n  1  n+1 

then  its  centroid  is  given  by 

n+1 

(22)  =  X  If"  • 

j  =  l 

Let  V  be  the  hypervolume  of  S  .  We  employ  the  superscript  n  in  our 
n  n 

notation  since  we  later  want  to  discuss  and  ^  simultaneously. 


THEOREM.  An  integration  formula  exact  for  the  general  cubic  polynomial 
over  for  n  >  1  is  given  by 

n+1 

(23)  f  f(x")  dx“  S  a  Y  £(x")  +  b  f(c“) 

n  j=l 

where 


(24) 


a  = 
n 


±21 


b  = 


V 


n+2)  n  ’  n  “  4(n+2)  n 


„n  2  „n  n+1  ^n  .  ,  ,  ^  , 

Xj  =  JiTj  ^  >  J  =  1  to  n  +  1  . 


(Before  beginning  the  proof  we  remark  that  the  points  are  on  the  median 

lines  of  S  and  the  statement  of  the  theorem  is  in  symmetric  form), 
n 


PROOF.  There  exists  an  affine  transformation  mapping  any  simplex  onto  any 
other.  Hence  to  prove  the  theorem  we  may  use  the  vertices 

y”  -  (0 . 0)  ,  Y^  =  (1,0,.,., 0) 

1  2 

. 0) .  . 

It  is  readily  verified  that  the  formula  given  holds  for  n  ■  1  and  n  ■  2. 

tl"*  1 

Hence  we  assume  that  It  holds  for  E  and  proceed  to  show  that  It  also 
holds  for  E*',  where  n  -  1  >  2.  Let  f (x^,x^, , , ,  ,x”)  ■  f(X*')  be  a  cubic 
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polynomial  In  x^,x^, . , . ,x".  Then  f ( l,x^, . . . ,x”)  is  a  cubic  polynomial 
2  n 

in  X  ,...,x  .  Now  using  a  result  established  for  cones  (equation  (I7))  we 


may  write 

(25) 


r  f(x")  dx”  =  rV  r  f(x")dx’^"^^  dx^ 

ds  do  \  J  I  J 

n  X  S  , 


’ (  l^^-l 

(x  ) 


n-1 

n+1 

n-lZ  Vl 

j=2 


dx 


n  —  1 

where  S  .  is  the  (n-l)  -  simplex  with  vertices  Y.  ,  Y 
n“  i  2 


n+1 


n-5-1 

(Y^  ^  =  Y^)  in  the  hyperplane  x^  =  1,  ^  ”  n  ^  where 

j=2 


X 


n-JL 

j  ’ 


^  and  b^  ^  are  given  by  equation  (24),  with  n  replaced  by  n-l, 

j  running  from  2  to  n+1.  It  is  observed  that  the  hypervolume  V  ,  is 

1  n-l, 

n  1  ^1  P 

l/(n-l)I  and  V  is  l/nl  .  Let  f(X  )  be  the  monomial  (x  )  (x*^) 

(x"^)  d  where  0  <  k,  +  +  k_  <  3.  Using  (25)  we  find  on  substitution 

=  1  2  3  = 

and  simplification  that 

2-k  -k  k  k  3"ko"k, 

r  1  h  n  ^(3  +3  ^+n“2)-n  ^  5; 

(26)  /  (x^)  ^  (x^)  2  (k3)  3  dx’"  =  -S.— 


J 


n 


4(n+l)  (n  +  k  +  k  +  k  ) 
123 


On  the  basis  of  our  assumption  (26)  gives  the  value  of  the  integral  indic¬ 
ated.  On  the  other  hand,  formula  (23)  (and  (24))  applied  to  f(X^)  = 

T  k  k  k 

(x)  (x)  (k^)  d  gives 

V  r  2-k  -k  -k  k  k  k  k  3-k  -k  -k  k  - 

(27)  4(;;;:nT^  ^  ^  ^[n  ^+(n+2)  ^^(3  ^+3  ^+n-2)]-(n+l)  ^ 


Now  it  may  then  be  directly  verified  that  (27)  gives  the  same 
result  as  (26)  for  0  <  ^  ^  ^  3>  ^2  -  ^3  ° 
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Hence  in  view  of  the  symmetry  with  which  the  last  n  ~  1  coordinates 

appear  in  the  set  of  vertices  Y  , . . . ,Y  (4)  is  verified  as  correct  for 

k/  ,  k 

all  monomials  of  the  form  (x  )  (x  )  -^  (x  )  for  0  <  k  -fk  +k,  <  3, 

-  1  J  *  = 

j  ^  /  where  j  and  I  are  taken  from  2  to  n.  The  only  monomial  type 

2  3  4 

thus  omitted  is  x  x"'^  x  provided  n  >  4o  Using  formula  (25)  for  this 
monomial  we  find 


n 


(28)  dx"  = 

n 

tiv  2  3  4 

which  coincides  with  the  value  obtained  on  substituting  f(X  )  =  x  x  x 
in  (25)  (and  (24)).  Hence  our  induction  is  complete  and  formulas  (23)  and 
(24)  hold.  (The  result  (23)  and  (24)  and  proof  is  due  to  P.  C.  Hammer  and 
A.  H.  Stroud,  published  in  M.T.A.C.  V.IO,  p. 


This  is  then,  a  useful  result  for  polyhedrons  in  that  every  poly¬ 
hedron  is  made  up  of  a  finite  number  of  n-simplexes.  However,  since,  for 
example  the  n-cube,  is  made  up  of  n!  n-simplexes,  it  is,  in  general,  worth¬ 
while  using  a  different  integration  formula  for  symmetric  polyhedrons. 


vi)  Thacher’s  matrix  method  of  attack 


It  was  noted  earlier  that  in  order  to  find  an  integration  formula 

of  the  form  (4)  of  precision  p  over  an  arbitrary  region  in  we 

would  need  to  solve  (n+l)p/p!  simultaneous  non-linear  algebraic  equations. 

For  w(X,)  !=  1,  these  equations  may  be  written  in  the  explicit  form 

,  m 

o  r  A  J  J 

(29)  I. 


S'  r,..  r  {nh''*- 

k^,  ...,kn  J 


isl 


dx^'  =  )  C  I, 
j=i 


ft 


n 


where  the  m  c.’s  and  the  mn  x/s  are  unknowns,  0  <  )  k,  <  p.  For  a 

J  j  “  i  “ 

1=1 


symmetric  region  the  number  of  equations  will  of  course,  be  considerably 
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1.  i 

less  since  we  need  make  no  distinction  between  the  variables  x  ,  x'^ 
(i, j=l, . , . ,n)  of  a  symmetric  region^ 


Define  a  set  of  m  x  m  diagonal  matrixes  by 


(30)  G  =  [G.  6^.] 


(51)  =  [xj  8^.1 


(52) 


In  terms  of  these,  matrixes  (29)  becomes 

ft 


tr  G 


k. 

1 


i=l 


=  I 


k.iy««»)k 

1  n 


In  particular,  for  hypercubes  of  edge  length  2. 


k, , . .  .  ,k 
1  n 


0  if  at  least  one  k,  is  odd 

1 

n 

n  (-—"■)  if  no  k.  is  odd. 

i=sl  Kj^+i  j 


For  a  second  degree  formula  we  have 


(55) 


tr  (G  G}  =  2” 


(5‘^) 

(55) 


tr  (GY^G)  =  0 


tr 


3  ij 


These  equations  in  the  traces  may  be  converted  to  vector  equations 

if  we  Introduce  the  m-dlmenslonal  column  vector,  e,  with  all  elements  unity, 

T 

and  its  transpose  e  .  Then,  if  we  define  the  vectors 

by  S  »  Ge,  q  «  Y^Ge.  «  Y^^Ge,  and  so  on,  (33),  (34)  and  (35) 

become 


G  G  s  .  (G  e)’^  (Ge)  -  -  2" 


(36) 


20 


(37)  G  G  €  =  (G  €)^  (Y^  g  e)  =  =  0 

(38)  G  Y^  Y-^  G  €  =  (Y^G€)^(Y^Ge)  =  =  (YS:^Ge)^(Ge)  =  =  |  6^. 

These,  however,  are  recognized  merely  as  orthogonality  relations 
among  (n-f-l)  vectors  ^^,..0,  and  the  normalization  requirement  that 

Kl^  =  2^  UJ2  ^  273.  Now  (n+1)  orthogonal  vectors  span  a  vector  space 
of  dimension  (n+l)  and  this  space  must  be  a  subspace  of  the  vector  space 
of  dimension  m  consisting  of  all  m-dimensional  vectors.  Thus  a  second 
degree  integration  formula  can  be  obtained  with  m  =  n+l  and  for  any  higher 
value  of  m. 


The  above  argument  has  also  furnished  an  explicit  algorithm  for 
constructing  examples  of  such  formulas  by  orthogonalizing  any  linearly 
independent  set  of  (n+l) (n+l)  dimensional  vectors,  and  applying  proper 
normalization  conditions. 


For  example,  orthogonalizing  the  set 
n/3  tan  6),  (3s/3  tan  9,  0,  O)  results  in 


(1,1,1).  (3,-  '/3  tan  0, 


(39) 

+0) 

(41) 


£  =  (2//3,  2//3,  Z/S) 

=  (™  cos  0,  ^  cos(0  + 

fc  /  2/2  .  _  2s/2  ,  /- 
^  =  (-—  sin  0,  sin(0  + 

23  3 


2jrv  ^ 

3^’  3 

y'  3 


cos(0  + 
sin(0  + 


—)) 

^)  » 
y* 


Hence,  for  this  case,  we  obtain  the  integration  formula 

...  r/  \  j  j  'v-  4  ft/™  cos  0,  sin  0) 

(42)  j  I  f(x,y)  dx  dy  ™  jf(^^ 

'-1  '-1 

+  f(^  cos(0  +  -“),  sln(0  +  ~)  +  f(^  cos(0  +  “), 
-P  sin-(0  +  “))j" 
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which  is  exact  whenever  f(x,y)  is  a  polynomial  of  degree  two  or  lees. 

For  3rd  degree  formulas  we  have,  along  with  (36),  (37)  and  (38) 
the  condition 


(1*5) 


G  yi  yJ  Y*"  G  e  =  (Y^  Y^  G  cf  (y''  G  c)  =  =  0 


Thus  we  must  consider  the  n{n+l)/2  new  vectors  in  addition  to  ^ 

and  These  new  vectors  fall  into  two  classes:  the  ortho¬ 
gonal  to  every  but  not  to  while  the  n(n-l)/2  (i  ^  j)  are 

orthogonal  to  both  sets,  or  else  are  null  vectors. 


Let  us  consider  first  the  case  where  all  the  t are  null  vectors.  This 

implies  that  unless  one  or  more  elements  of  ^  is  zero,  in  which  case  the 

basic  integration  formula  includes  redundant  points  with  zero  weight,  only 

one  of  the  can  have  any  given  element  different  from  zero.  The 

cannot  be  null  vectors  in  view  of  (38)  while  from  (37)  unless  the  differ 

± 

in  sign,  Y  must  include  elements  of  both  signs,  and  thus  at  least  two 
non-zero  elements.  We  therefore  conclude  that  the  dimension  m  of  Y^  must 
be  at  least  2n. 


For  an  equally  weighted  formula  G  is  a  scalar  matrix  which  we 

may  write  as  gE  where  E  is  the  identity  matrix.  From  (33)  fo*'  ®  2n 

n*  1 

point  formula,  g®  must  have  the  value  2  /n.  For  the  minimum  number 

of  points  Y  will  have  but  2  non-zero  elements  of  opposite  sign  and  from 
(34)  of  equal  magnitude,  and  these  may  be  arranged  in  order  such  that 

(kk)  . 

a  diagonal  matrix  with  the  (2i-l)’th  element  equal  to  x^  and  the  (2l)'th 
to  -x^,  From  (35)  it  follows  that 
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(45)  2(x^)2  g2  =  2n/3 

so  that  for  all  i 

(1.6)  ^ 

and  we  have  the  family  of  2n  point  third  degree  integration  formulas 
given  in  equation  (lO). 

If  one  or  more  of  the  t,.  is  non-null  there  must  be  n+2  or 

ij 

more  orthogonal  vectors,  and  the  minimum  value  of  m  must  be  at  least  n+2. 
This  lower  bound  cannot  always  be  achieved,  the  only  region  for  which  we 
have  seen  an  (n+2)-point  third  degree  formula  being  the  n-simplex  (formula 
(23)).  Thacher  claims,  moreover,  that  it  can  be  shown  that  no  5  point 
third  degree  formula  exists  for  the  cube. 

This  method  of  obtaining  integration  formulas,  due  to  H.  Thacher 
Jr.  is  interesting  in  that  he  outlines  a  different  attack  by  which  he  first 
finds  the  minimum  number  of  points  required  and  then  proceeds  to  find  the 
points  and  weights  (V.ll  p.l89  M.T.A.C.).  For  higher  precisions  than  3  the 
method  becomes  exceedingly  complicated,  since  additional  non-linearities 
are  introduced. 


vii)  A  remark  on  the  uniqueness  of  evaluation  points 

Note  that  for  n  >  ^  in  formula  (lO)  of  Chapter  I,  the  points 

x^  lie  outside  the  n-cube.  Since  the  distance  from  the  center  of  the  n- 
J 

cube  to  its  furthest  edge  is  ^yn,  sfn>  ,  so  that  for  n  >  3  we  can 
rotate  the  points  of  Integration  back  into  the  n-cube,  provided  a  rotation 
of  the  points  does  not  reduce  the  precision  of  our  integration  formula. 


Now,  a  rotation  is  a  particular  affine  transformation  for  which 


the  elements  of  the  matrix  the  transformation  satisfy 
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n 


(4?)  y 


i=l 


When  integrating  the  monomial  x^  of  degree  2  over  the  n-cube,  we 


obtain 


r'-  -k  -t  -1  -n  a"  . 

I  ...  I  X  X  dx  . .  .  dx  =  —  o,  . 

J  J  3 

-1  -1 

n 


-k  i 

Applying  the  rotation  x  =  y  x  and  integrating  the  transformed 


i=l 

monomial  over  the  same  n-cube  we  obtain 

n  n 


(^9) 


-1  -1  i=l  j=l 

n 


=  f\..  r 

J  J 


1  r 


-1  -1 


a 


ik 


i=l 


(I  “ik  “j/ 

3 


Hence  we  have  shown  another  way  of  justifying  a  result  for  the 
n-cube  due  to  A.  H.  Stroud  (V.ll  p.257»  M.T.A.C.),  if  Xj  a  (Xj,,,.,Xj) 

21in 


then 

(50) 

21-1 

1  *?'. 

9  008  a 

3  n4'l 

g  Bin 
3  n+1 


(If  n  is  odd,  Xj  ■  i~l)^/J3) 


Moreover  we  have  shown  that  for  the  n-cube  the  points  of  inte¬ 
gration  of  the  second  and  third  degree  integration  formulas  are  not  uni¬ 
quely  determined. 
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viii)  A  precision  2  Tonnula  over  an  arblLrary  region  lii 

A.  H.  Stroud  (V.l4,  p.21,  M.T.A.C.)  extended  H.  Thacher's  matrix 
attack  to  obtain  an  n+1  -  point  formula  of  precision  2  valid  over  an  arb¬ 
itrary  region  in  e”  that  (together  with  the  arbitrary  weight  function) 
satisfies  certain  conditions  of  non- degeneracy .  We  summarize  his  proced¬ 
ure  in  what  follows. 


Assume  the  integration  formula 

n 

w(X)  f(X)  dX  ^  y  c.  fix.) 
n  j^O 

to  be  of  precision  2  and  to  require  n+1  points, 

at  this  point  we  define 


To  establish  our  notation 


(51) 


a  =  r  w(X)  dX 

oo  Jr  '  ' 

n 


a  =  r  w(X)  x^  dX 
oi  Jr 

n 


a,  ,  =  r  w(X)  x^  dX 
Jr 

n 


To  obtain  the  above  integration  formula  we  want  to  find  the 

numbers  c,  and  the  points  v^  given  in  the  equations 
J  J 

n 

V" 

)  c  .  =  a 

L  J  oo 

j=0 
n 


(52)  < 


V  1 

)  c  .  V ,  =  a  , 
A  J  J 


j=0 


n 


V’  k  i 

>  c  ,  V  ,  V .  =  a 

U  J 
j=o 


j  j 


We  can  write  (52)  in  the  matrix  form 


(55)  V  C  V  =  A 


25 


where 


(51*)  \ 


V  = 


n 


V, 


n 


'^0 

0 

0* 

00 

^01 

^On 

c  = 

0 

"l 

0 

,  A  = 

01 

a^l  ... 

^In 

0 

0 

*  •  •  c 

n 

On 

^In 

a 

nn 

and  where  we  assume  0  <  a^^  <  oo  and  A  non-singular. 

00 

Since  A  is  symmetric  we  can  always  find  a  non-singular  transfor¬ 
mation  B  such  that 


(55)  C  V  B  =  A  B  =  a^^  D 


where  D  is  the  diagonal  matrix  with  elements  +  1. 
We  set  VB  =  Z  ,  where 
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Post -mult ip lying  (55)  by  D  ^(VB)^  we  obtain 
(VB)^  C  (VB)  D"’-  (VB)'^  .  (VB)^ 


(57) 


I 

I 


\  V. 


(r?> 


a 
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Pre-multlplying  (5?)  ^  [(VB)  ]  (C  is  non-singular)  and  noting  that 

D  ^  =  D  we  obtain 


(58)  (VB)  D  (VB)*^  =  a^Q  c"^  =  Z  D  . 


In  terms  of  the  this  equation  is 


(59) 


'  ^  i  \  -  I  \  = 


i=l 


i=r+l 


00  5- 

c  ki 
k 


where  r+1,  0<r<n  is  the  number  of  +  I's  in  D. 


We  are  only  interested  in  real  solutions  of  (52)  and  therefore 
precisely  n  -  r  +  1  of  the  must  be  negative  by  Sylvester's  "law  of 

inertia"  ([9]  P-56).  Keeping  this  in  mind,  we  first  find  a  particular  set  of 
z^  's  that  will  satisfy  equations  (59)-  We  then  use  equation  (56)  to  find 
the  inverting  the  matrix  B;  i.e. 


(60) 


V  =  Z  B 


-1 


From  a  particular  solution,  other  solutions  may  be  found  as  follows 


If 


(61) 


S  = 


1  0 

0  s 

0  s 


11 


0 

s 


nl 


in 


nn 


T 

is  a  cogredient  automorph  of  D,  that  is  if  SDS  =  D  then  by  equation  (58) 
we  see  that 


(62) 


Y  =  Z  S 


can  be  used  Instead  of  Z  in  equation  (60)  to  obtain  another  solution  for  V. 
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Table  I  gives  a  particular  solution  of  (58);  we  have  assumed 
^0*’’’’*^n-r  and  ^n+l-r  ’  ’  ‘  ’  ^n  Positive.  In  cases  where  the  double 

sign  occurs  we  mean  to  use  the  upper  sign  for  the  first  r  components  of 

1  ri 

each  vector  =  {2.y...,z^)  and  the  lower  sign  for  the  last  n  -  r  com¬ 
ponents.  Each  Zj  is  real. 

TABLE  I 
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The  formulas  discussed  above  are  minimal;  that  is,  similar  formulas 
cannot  be  obtained  with  m  +  1  points,  where  m  <  n.  For  if  a  formula  could 


Mil-' 
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be  obtained  with  m  +  1  points  ,  j=0,...,ni,  m  <  n  then  equation  (53) 

would  still  hold  where  A  is  the  same  as  before  and  U  is  (m+l)x(n+l),  C 

T 

is  (m+l )x(m+l ) .  Hence  det  (U  CU)  =  0  since  det  U  =  0.  By  assumption 
det  A  ^  0  and  thus  (53)  cannot  hold  if  m  <  n. 


Proceeding  in  this  manner,  A.  H.  Stroud  also  obtained  a  2n-point 
precision  3  integration  formula  valid  over  a  symmetric  but  otherwise  arbit¬ 
rary  region  in  In  this  case  the  weight  function  is  also  symmetric,  but 

otherwise  arbitrary.  In  the  same  paper  he  gAy^  results  for  the  precision  5 
case  similar  to  those  in  Table  I  and  he  also  gAv'e  methods  for  extending  these 
results  to  other  regions  in  e”  by  linear  transformations. 


ix)  Monte  Carlo  method 

A  very  well  known  method  of  numerical  integration  applicable  over 
an  arbitrary  region  in  e"  is  the  Monte  Carlo  method.  This  method  essenti¬ 
ally  consists  of  the  evaluation  of 

m 


(63) 


r  f(x)  dx  =  -  y  £{x.)  ;  V  =  r 

Jr  '  '  ”>  L  y  Jr 

n  j=l  n 


dX 


where  the  X^'s  are  points  chosen  uniformly  at  "random”  throughout  the 
region  of  integration,  R^, 


If  f(X)  is  sufficiently  well  behaved  (e.g.  of  bounded  variation 
in  the  variables  X)  it  can  be  shown  statistically  that  the  error  of  inte¬ 
gration 


(6M 


We  have  quotations  around  the  word  random  above  since  the  points  Xj  are 
obtained  by  the  use  of  a  formula,  for  example 
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^i+1  ^ 

As  an  example  of  the  accuracy  obtainable  by  this  method,  Davis  and  Rabinowitz 
(Vol.  10,  p.l,  M.T.A.C.)  show  they  can  obtain  results  to  within  one  per  cent 
error  using  approximately  10,000  points  to  compute  the  volume  of  a  4-dimen¬ 
sional  sphere.  Other  examples  are  given  on  page  48  of  this  thesis. 

Suppose  that  due  to  a  large  value  of  n  and  complexity  of  f(X) 
it  takes  an  I.B.M.  7^90  one  hour  to  compute  a  solution  accurate  to  one  sign¬ 
ificant  figure.  To  achieve  4-figure  accuracy  would  then  take  over  20  years. 

Due  to  convergence  to  an  accurate  solution  being  slow  by  the  Monte 
Carlo  method,  the  method  is  used  normally  to  make  rough  estimates  of  integrals 
of  a  dimension  so  high  that  it  would  be  impractical  using  a  numerical  inte¬ 
gration  formula. 

For  a  detailed  explanation  of  Monte  Carlo  methods  see  the  paper 
"Monte  Carlo  methods  for  solving  multivariable  problems",  by  J.  M.  Hammers- 
ley  in  Annals  of  the  N.Y,  Acad.  Sc.  I96O  Vol.  86,  Art.  3,  p.  845”874.  This 
paper  also  contains  a  wealth  of  references  on  the  method, 

x)  Taylor  Series  approach 

J.C.P.  Miller  has  written  several  papers  ( [lO] , [ 11 ] , [ 12] )  on  inte¬ 
gration  over  a  rectangular  domain.  To  obtain  an  integration  formula  of  the 
form 

h  h  ™ 

(65)  1=  f  ...  f  f(X)  dX  S'  (2h)"  )  c  £(X  ) 

J  J  ■}  J 

-h  -h  2=1 

he  expands  an  arbitrary  function  f(X)  in  a  Taylor  series  in  n  dimensions 
and  tries  to  chose  his  m  points  over  the  region  of  Integration  (the  n-cube) 
in  such  a  way  as  to  represent  the  function  by  as  many  of  the  first  terms  of 
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the  Taylor  series  expiinsion  as  possible. 


Integrating  his  Taylor  series  expansion  over  the  n-cube  of  side  h 


he  obtains 


(66)  I  -  (2h)"  |£q+j,  Sg+j,  (S,^+S2_2)+^;  *®8'^®6,2 

"^^4  4*^^ 4  2  2*^^ 2  2  2  2^"^  *  * 

where  by  we  mean  f(0,...,0)  and  where,  for  example 


(67) 


S  =  s 

2,2,2  2 


A* 

.2.2  L  7D- 


f. 


i,/rk=i 


the  asterisk  indicating  that  i,  j,  k  are  unequal  in  pairs  throughout. 


®2,2.2  “  5'  5'  • 


Setting  up  and  solving  non-linear  algebraic  equations  he  obtains 
various  integration  formulas;  some  of  these  being  similar  to  those  that  have 
previously  been  obtained  by  Hammer  and  Stroud,  and  which  we  will  list  later. 
His  procedure  also  gives  us  an  error  estimate  depending  on  high  order  partial 
derivatives  of  the  function  f(X),  and  although  such  an  error  estimate  may 
be  difficult  to  evaluate,  his  papers  are  among  a  very  few  that  discuss  errors 
of  numerical  Integration  formulas  in  higher  dimensions  at  all. 

Moreover,  by  his  method  he  is  able  to  see  that  harmonic  functions 
require  fewer  points  of  evaluation  than  arbitrary  functions,  and  he  exploits 
this  in  [12]  using  polynomials  orthogonal  over  the  Interval  (0,l)  with 
reapect  to  weight  functions  -  x  -  x^^^)  .  As  an 

example  he  evaluates 
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to  an  accuracy  of  ten  significant  figures  using  only  five  points  on  the  square 
(We  recall  here,  that  both  cos  x  and  cosh  x  have  very  riipidly  converging 
Taylor  series  expansions,  and  hence  this  is  probably  an  example  of  the  Integra 
tion  formula  at  its  best.  The  accuracy  of  the  result  is  surprising,  neverthe¬ 
less)  . 


xi)  The  error  bound  of  von  Mises 

We  close  this  chapter  with  a  summary  of  an  error  bound  due  to  von 
Mises  [26],  and  summarized  by  Stroud  [27]. 


Suppose  we  have  the  integration  formula 


m 


(68)  r  £(X)  dX  =  y  c.  f(X.)  +  E 

Jr  £_j  J  J 

n  j=l 


where  R  is  starlike*  with  respect  to  an  interior  point  Z, 
n 

to  spherical  coordinates,  von  Mises  shows  that 

rn 

L_  r  P  r*^  dX  +  y  jc,]  r^ 

• '  I  I  i  R  ^ 


Transforming 


(69) 


El  <  If 


'max  (i; 


j  =  l 


where 


r5=  [(x^-z^)^  +  ...  + 


f 


(M)  _ 


dx 


1 


r  n 
oz  J 


f(X) 


Here  if  the  integration  formula  has  precision  p  then  p  can  be 
taken  to  be  l,2,...,p+l  ;  r  is  the  distance  from  Z,  and  r^  is  the  dis¬ 
tance  of  the  j'th  point  from  Z. 


*  A  region  in  I0  star  like  with  respect  to  an  interior  point 

Z  if  any  half  ray  emanating  from  Z  cuts  the  boundary  of  at 

one  and  only  one  point. 
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For  even  values  of  n  <  p  the  Integral  /  r^  dX  will  be  known 

n 

from  the  monomial  integrals  required  to  determine  the  integration  formula. 
For  odd  values  of  this  integral  may  be  difficult  to  evaluate,  and  in 
such  cases  it  may  be  desirable  to  use  the  Schwarz  inequality 


where  V(R  )  is  the  volume  of  R 
n  n 


If  the  c  are  all  positive,  then  von  Mises  shows  that  equation 

j 

(69)  can  be  improved  to 


(71) 


m 


f  r-dX. 

.  f(^)  y 

L  V. 

min  / 

n 

m 

r  r*^  dX  - 
'  R 

.  y 

L  ^ 

max  /  , 

n 

Equations  (69)  and  (7I)  are  valuable  theoretical  results,  but  find¬ 
ing  the  derivative  of  a  complicated  integrand  and  obtaining  a  bound  on  it  can 
be  extremely  difficult  in  practical  cases,  particularly  for  large  n  and/or 

P. 


xll)  Summary 

The  theory  and  results  presented  up  to  this  point  constitute  a 
rough  outline  of  current  methods  of  attack. 


To  restate  the  problem,  we  should  like  to  have  numerical  integra¬ 
tion  formulas  of  the  type  (1|.)  of  Chapter  I  which  employ  values  of  the  inte¬ 
grand  at  points  lying  within  the  region  of  integration.  The  formulas  are  to 
have  precision  p  as  described  on  page  3  total  number  of  points  m 

is  to  be  as  small  as  possible. 
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The  solution  of  this  problem  faces  us  with  another  problem;  the 
solution  of  a  set  of  simultaneous  non-linear  algebraic  equations*.  Hence 
if  effective  methods  of  solving  a  set  of  simultaneous  non-linear  algebraic 
equations  were  available,  the  task  of  finding  numerical  integration  formulas 
in  higher  dimensions  would  not  be  such  a  difficult  problem. 

In  setting  up  these  non-linear  algebraic  equations  it  would  help 
considerably  if  we  knew  in  advance  the  minimum  number  of  points  m  =  m(n,p,R^,w) 
required  to  obtain  integration  formulas  of  the  type  (h)  of  Chapter  I.  We 
already  know  from  a  study  of  integration  formulas  in  one  dimension  that  this 
minimum  number  of  points  may  be  dependent  on  w(X).  There  is  some  evidence 
that  this  minimum  number  of  points  may  also  be  region  dependent.  To  illust¬ 
rate  this  point  we  have  presented  above  the  derivation  of  an  n  +  2  -  point 
precision  5  integration  formula  for  the  n-;&implex  (Hammer  and  Stroud  [6]), 
while  the  best  existing  formula  of  precision  5  uses  2n  points  (it  is  not 
known  whether  this  is  a  minimal  formula).  Thacher  [22]  claims,  moreover, 
that  a  5"Point  precision  5  formula  for  the  cube  does  not  exist. 

It  is  desirable  to  keep  n  arbitrary,  since  we  wish  to  achieve, 
general,  rather  than  particular  integration  formulas  in  higher  dimensions. 
Although  we  lose  somewhat  In  generality**  if  we  restrict  our  considerations 
to  fully  symmetric  regions  we  gain  in  algebraic  convenience  since  we  shall 
have  fewer  non-linear  equations  to  solve  and  the  number  of  these  equations 
will  not  depend  on  n. 


*  Workers  in  research  are  faced  with  the  non-linear  problem  far 

more  frequently  than  they  are  faced  with  the  linear  one. 


** 


The  integration  formulas  may  undergo  a  loss  of  precision  under 
transformation  (transformation  theorem,  Hammer  et  alii  p.9)' 
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Error  estimates  are  necessary.  In  one  dimension,  the  error  of 
integration  formulas  is  popularly  expressed  in  the  form  of  a  high  order 
derivative  of  the  integrand--obtaining  a  bound  on  an  error  estimate  such 
<18  this  may  be  extremely  difficult,  and  this  is  much  more  so  in  higher 
dimensions . 
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CHAPTER  III 

THE  CONSTRUCTION  OF  INTEGRATION  FORMULAS 

In  this  chapter  we  compare  the  usefulness  of  formulas  which  inte- 

n  .  k. 

integrate  all  monomials  of  the  form  IT  (x  )  ,  k,  <  p  (the  formulas  that 

i=l  ^  ^  “ 

are  Cartesian  products  of  formulas  in  one  variable).  The  minimum  number  of 
points  required  to  obtain  an  integration  formula  of  precision  p  over  an 
arbitirary  region  in  n-space  is  estimated. 

Starting  with  appropriate  sets  of  symmetrically  spaced  points  non¬ 
linear  algebraic  equations  are  set  up  and  solved  to  obtain  fully  symmetric 
numerical  integration  formulas. 

Orthogonal  polynomials  are  constructed,  and  as  for  the  case  of  one¬ 
dimensional  integration  the  zeros  of  these  polynomials  (at  least  the  zeros  of 
the  odd  polynomials!)  turn  out  to  be  the  required  non-linear  unknowns  in  the 
simultaneous  non-linear  algebraic  equations  set  up  to  obtain  numerical  inte¬ 
gration  formulas.  The  polynomials  were  thus  an  aid  in  obtaining  the  solution 
to  the  non-linear  algebraic  equations. 

l)  A  Discussion  of  Precision 

We  recall  (see  Cartesian  product  theorem  Chapter  II,  p. ll)  that 
high  precision  integration  formulas  are  available  over  Cartesian  product 
regions  in  n-space.  For  example,  to  obtain  a  formula  of  precision  2m-l  we 
would  need  m  points  on  a  line,  m^  points  on  a  rectangle,  ...,  m”  points 
on  a  rectangular  region  in  n-space.  To  obtain  a  formula  of  precision  5  over 
the  n-cube  we  would  need  2^^  points,  while  formula  (lo)  of  Chapter  I  tells 
us  that  we  can  obtain  precision  3  with  only  2n  points. 

This  is  a  somewhat  unfair  comparison  however,  since  a  repeated 
Gaussian  integration  formula  will  integrate  so  many  more  monomials  than  an 
Integration  formula  which  will  only  integrate  monomials  up  to  degree  p.  A 
repeated  Gaussian  integration  formula  will  Integrate  any  monomial  of  the  form 


grate  all  monomials  of  the  form  |~[  (x  ) 

i=l 


y 


^i  ^ 


p  with  formulas  which 
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(l)  H  (x^)  ^  where  0  <  <  P 

i=l  -  1  - 

This  monomial  is  of  degree  <  p  in  each  variable  and  of  total  degree  <  n  p  . 


The  definition  of  precision  (an  integration  formula  is  said  to  have 

fr  i  ^i 

precision  p  if  it  exactly  integrates  any  monomial  of  the  form  (x  ) 

n  i=l 

where  0  <  ^  ^  p)>  then  has  made  a  whole  series  of  integration  formulas 

iTl 


possible;  formulas  which  at  first  notice  seem  to  perform  as  well  as  the  re¬ 
peated  Gaussian  formulas,  but  require  by  far  fewer  points.  Now  it  is  true 
that  the  formulas  arising  as  a  result  of  this  definition  of  precision  require 
fewer  points  than  the  repeated  Gaussian  integration  formulas,  but  this  is 
mainly  because  the  former  formulas  integrate  fewer  monomials  than  the  latter. 
With  a  precision  3>  2^  "  point  repeated  Gaussian  integration  formula  we 
can  exactly  integrate  U”  different  monomials,  while  with  the  2n  -  point 
formula  (lO)  of  Chapter  I  we  can  exactly  integrate  (n+l) (n+2)(n+3)/3-’  diff¬ 
erent  monomials.  For  n  =  20,  the  ratio  of  the  number  of  monomials  exactly 

12  3 

integrated  by  the  two  formulas  is  roughly  10  to  2x  10  . 


It  is  also  interesting  to  compare  the  number  of  monomials  each  of 
these  formulas  integrate  per  point.  Any  repeated  Gaussian  integration  form¬ 
ula  will  integrate  2^^  monomials  per  point,  while  formula  (lO)  of  Chapter  I 
will  integrate  only  0(n  )  monomials  per  point.  We  conjecture  that  2 
monomials  per  point  is  the  best  possible. 


(2) 


The  polynomial 
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-  ^ 

■j'^0 
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0  <  k,  <  P 

=!  1  =S 


has  monomials  up  to  degree  p  in  each  independent  variable.  That  is,  it 
contains  all  the  monomials  an  arbitrary  "Cartesian  product"  type  polynomial 
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of  degree 

(3) 


p  would  contain.  The  polynomial 


P  JL 


n 

[a,  .  V 

i 

a ,  X 

j=i 

0 

u 

it 

h  J 

contains  all  the  monomials  an  arbitrary  Taylor  polynomial  of  degree  p  would 
have . 


Now  any  linear  or  affine  transformation  will  not  alter  the  degree 
of  (3)«  However,  a  linear  or  affine  transformation  will  increase  the  degree 
of  (2),  unless  the  transformation  is  of  the  form  =  a^x^  +  .  This 

seems  to  indicate  that  the  integration  points  for  a  repeated  Gaussian  inte¬ 
gration  formula  are  uniquely  determined.  It  suggests,  moreover  (as  has  al¬ 
ready  been  shown  for  a  particular  case  on  page  23)  that  the  points  x^  of 
integration  formulas  capable  of  integrating  only  the  monomials  up  to  degree 
p  (and  not  products  of  monomials  up  to  degree  p  in  each  independent  vari¬ 
able)  are  not  uniquely  determined. 


When  integrating  a  monomial  of  the  type  (5)  of  Chapter  I,  we  inte- 

i 

grate  with  respect  to  each  variable  x  holding  the  remaining  variables 
fixed.  Hence  in  this  application  we  are  concerned  more  with  the  highest 
degree  in  each  variable  than  with  the  total  degree,  in  all  the  variables,  i.e. 


the  degree 


k. 

1 


in  each 


i 

X 


is  more  significant  than  the  total  degree 


i=l 


All  the  above  arguments,  together  with  the  difficulty  of  obtaining 
interpolation  polynomials  that  represent  an  arbitrary  polynomial  of  degree 
p  In  n  variables  (see  H.  C,  Thacher  Jr.,  New  York  Acad.  Sci.,  v.86, 
art. 3,  p. 758-775)  seen  to  indicate  that  the  definition  of  precision  given  on 
page  3  of  this  thesis  does  not  fit  too  well,  particularly  for  integration 
formulas  over  rectangular  regions  in  hyperspace. 
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We  must,  nevertheless,  not  forget  that  the  total  number  of  points 
required  in  a  repeated  Gaussian  integration  formula  becomes  fantastically 
large  as  n  increases.  It  may,  moreover,  often  occur  that  we  can  obtain 
comparable  accuracy  using  a  formula  that  requires  by  far  fewer  points.  To 
illustrate  this,  let  us  estimate  the  number  of  points  m  required  to  obtain 
an  integration  formula  of  precision  p  over  an  arbitrary  region  in  n-space. 


We  have 
dinates  to  fit  an 

,<n+pv  .  . 

{  )  monomials. 

P 


at  our  disposal  a  choice  of  ra 
arbitrary  polynomial  of  degree 
Hence  m(n+l)  >  or 


m  > 


pl  (n+l)I 


weights  and  mn  coor- 
p  in  n-space  that  has 


When  n 

is  large  we  may  use  Stirling's  formula 

(5) 

r(n+p+l)  'x.  n^‘*'P+2  g 

to  obtain 

nP-^ 

(6) 

m  >  , 

=  p: 

as  an  estimate  of  the  minimum  number  of  points  m  required  to  obtain  an 
integration  formula  of  degree  p  over  a  region  in  n-space*.  The  repeated 
Gaussian  formulas  require  points. 


In  the  Table  on  the  following  page  we  have  tabulated  the  minimum 
number  of  points  required  by  a  formula  of  precision  p  along  with  the  number 
of  points  required  by  a  Gaussian  formula  of  precision  p. 


For  an  improvement  of  this  bound  to  ra  >  n  /kl ,  k  an  integer, 

P  =  2k  sec  [27];  the  fully  symmetric  formulas  of  the  following 

2k  ^^2k  ^ 

sections  require  m  =  points  for  large  n 


(2k) 


and  precisions  p  =  4k+l, 
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TABLE  II 


n 

p 

Gaussian 

riinimum  required  p  ( n  1 )  * 

k 

5 

26 

81 

4 

9 

143 

625 

4 

11 

273 

1,296 

4 

19 

1,771 

10,000 

8 

9 

2,601 

390,625 

8 

19 

247,000 

100,000,000 

This  Table  clearly  indicates  that  the  Gaussian  integration  formulas 
reach  the  economically  feasible  point  much  sooner  than  formulas  that  use 
a  minimum  number  of  points  for  even  the  fastest  computers.  We  would  there¬ 
fore  like  to  suggest  that  if  we  had  "minimum  point"  formulas  we  could  in 
most  cases  obtain  accurate  results  with  these  at  a  lesser  cost  than  with 
repeated  Gaussian  formulas. 


2)  The  Setting  up  and  Solution  of  Non-linear  Algebraic  Equations 

The  approach  given  here  is  similar  to  that  given  in  [7]  by  Hammer 
and  Stroud.  In  looking  for  fully  symmetric  integration  formulas*  we  first 
chose  appropriate  sets  of  symmetrically  spaced  points.  The  coordinate  sets 
and  the  number  of  points  obtainable  with  each  are  tabulated  below.  The 
method  of  counting  these  is  as  follows. 


Suppose  we  have  k  different  cocjrdinates :  r.  of  the  first, 

A  ^ 

of  the  kth  so  that  ^  a  n,  Then  by  rearranging  these  coordinates 
in  all  possible  ways  we  can  obtain  n!  / 


.5. 


different  points  in  n-space. 


If  we  further  allow  all  possible  changes  of  signs  of  the  coordinales  we  can 
k 

obtain 


j 


different  points  in  n-space,  where  w©  have  assumed 


For  the  definition  of  a  fully  symmetric  integration  formula  see 
page  8  . 


q  <  n  non- zero  coordinates. 


Coordinates 

Restriction  on  n 

Number  of  Points 

(0, , . . ,0) 

n 

> 

1 

2°(S) 

(+u,0, . . . ,0) 

n 

> 

1 

(±u,+u,0, . . , 

,0) 

n 

> 

2 

2^(2) 

(+u,±v,0, . . . 

,0) 

n 

> 

2 

2^(2)2! 

(+u,+u,+u,0, 

n 

> 

3 

23(5) 

(±Uj±t»,+U,+U 

0 

<D 

n 

> 

h 

2‘‘(!;) 

Stroud  in  v,l4,  p.21,  M.T.A.C.  gives  precision  5  integration  form¬ 
ulas  valid  over  an  arbitrary  symmetric  region  in  Note  the  distinction 

between  symmetric  and  fully  symmetric,  A  cylinder  for  example  is  symmetric, 
while  the  n-cube  is  both  symmetric  and  fully  S5mimetriCo 

The  integration  formulas  we  will  try  to  obtain  are  for  fully  symmetric 
regions.  The  weight  function  of  these  integration  formulas  is  also  assumed 
to  be  fully  symmetric,  but  otherwise  arbitrary. 

From  the  definition  of  fully  symmetric  it  follows  that: 

(a)  the  integral  over  a  fully  symmetric  region  R  of  any  product 
of  the  coordinate  variables  which  contains  an  odd  pox-rer  is 
zero;  and 

(b)  the  integral  of  a  product  of  even  powers  depends  only  on 
their  exponent  and  not  on  their  ordering. 

In  what  follows  we  have  written  down  non-linear  algebraic  equations 
and  their  solutions  to  obtain  fully  symmetric  integration  formulas  of  preci¬ 
sion  5,  5  and  9,  The  method  of  solving  the  non-linear  algebraic  equations 

is  given  in  the  next  section.  On  the  right  hand  side  of  the  non-linear  alge- 

riL  2  "i  ^ 

w(X)(x^)'^(x^)^  dX,  i  /  j  for  example  by 


41 


l2^(=  assume  also  that  contains  the  origin*, 

i)  Formulas  of  precision  5 

To  obtain  the  2n-point  equi-weighted  integration  formula  of  prec¬ 


ision  3> 


1(f) 


2n 

r  w(X)  f(X)  dX  =  A  y  f(+u,0, . . . ,0)  , 
Jr  Lj 

n  j  =  l 


we  solve  the  following  pair  of  simultaneous  non-linear  algebraic  equations; 


(7) 


f  £n  A 
I'  2A  u^  = 

2 


The  first  equation  here  is  obtained  by  substituting  the  integral 
of  a  constant  into  the  above  integration  formula;  the  second  by  substituting 
the  integral  of  the  square  of  a  particular  variable. 


General  Solution.  We  easily  obtain  the  value  of  A  from  the  first 
of  the  above  non-linear  equations;  the  second  equation  then  gives  us  the 
value  of  u.  The  general  solution  is 

Four  particular  solutions  are  tabulated  in  Table  III. 


il)  Formulas  of  precision  5 


We  assume  that  we  can  obtain  a  2n^  -i"  1  "  point  formula  of  the  form 
2n  2n(n“l) 

1(f)  »  A  f(0 . 0)  +  By  f(±u,0 . 0)  +  Cy  f(±u,tu,0,...,0) 

1”'  T 


Due  to  our  choice  of  a  symmetric  set  of  points  all  odd  powers  of  monomials 


* 


It  is  felt  that  this  assumption  may  not  be  necessary. 
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will  sum  to  zero.  To  obtain  the  values  for  u,  A,  B,  C 
following  system  of  non-linear  algebraic  equations 


(8)  < 


1 

2n 

2n(n-l) 

0 

2u^ 

4(n- 

•Du 

0 

2u^ 

4(n- 

.l)u 

0 

0 

4u 

we  need  to  solve  the 


General  Solution.  The  value  of  u  is  easily  obtained  by  dividing 
the  second  equation  into  the  third.  By  ordinary  elimination  we  can  then 
obtain  the  value  of  the  linear  coefficients.  The  general  solution  is 


A  - 

B  =  i  (yi^)2  [1^  -  („-l) 
°  =  \  ( 


Note  that  we  have  assumed  n  >  2,  Four  particular  solutions  are 
tabulated  along  with  the  precision  3  formulas  in  Table  III. 


iii)  Formulas  of  precision  9 

Due  to  difficulties  involved  which  will  be  made  clear  in  the 
following  sections  we  have  not  obtained  the  formulas  of  precision  7.  We 
have,  however,  obtained  the  formulas  of  precision  9«  the  presentation 

of  these  formulas  we  use  the  notations 

n^^^  n(n-l) .  .  ,  (n-k-i-l) 

^"(k)  "  k! 
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ji.i 


■jj 
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r(l+n/2) 


..  i 


n 


ta 

1  M 
}  M 

1  v-( 

H 

iX 

,..a 

1  ^ 

5 

A> 

i  C'<i 

>• 

y 

! 

1 

cf 

k! 


{t.  - 


^i  'ji  il  j' 


The  formulas  are  assumed  to  have  the  form 


2n 


1(f)  =  A  f(0,...,0)  +  y  [B  f(iu,0, . . . ,0)  +  C  f(±v,0,...,0)] 


2^n 


(2) 


y  [D  f(±u,±u,0, . . . ,0)  +  E  f (+v,+v,0, . . . ,0) ] 
+  y  F  f  (+u,+^v,0,  ...  ,0) 


2\ 


(J) 


y  [G  f (±u,+u,+u,0, . . . ,0)  +  H  f (+v,+v,+v,0, . . . ,0) ] 


2^n 


ih) 


y  [I  f (±Uj+u,+u,+u,0, . . . ,0)  +  J  f (+v,+v,+v,+v,0, . . , ,0) ]  . 


To  obtain  the  values  of  u,  v,  A,  . . . ,  J  of  this 
■~(n^  -  4n^  +  lln^  -  7n)  +  1  -  point*  formula  we  need  to  solve  the  follow¬ 
ing  system  of  non-linear  algebraic  equations: 

(9)  See  page  45* 


General  Solution.  This  is  valid  provided 


(10) 


Ig2  +  “^''^^22  =  ^1.-=  **  • 


h2 


The  method  by  which  this  solution  was  obtained  is  given  later. 


(11 


)  u,v  =  p 


2  8  4  6  X  L  v-^2-8^  "^^  4  2  6  ^^2"4"6"8 


^62  ~  ^44 

,2  2/  2  2v2 

4u  V  (u  -V  ; 


*G  =  0  This  reduces  the  number  of  integration  points  to 
^(n^'  -  5n^  +  I4n^  -  7n)  +  1  . 


** 


So  far  xve  have  found  no  examples  for  which  this  equation  does 
not  hold. 
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H  .  Ha2  ~  ^"'^^^2222 

■  8v8 

^422  ~  ^  ^222 

I6(n-3)u^(u^-v^) 


I  = 


J  = 


^■2222  " 


8 


l6v 


8 


^  '  ~'lt'f 2"  2**^  ■  -  2("-2)  [H  +  (n-3)J] 

4v  (u  -V  ) 

°  '  2  ^2v  ■  -  2(n-2)(n-3)l 

4u  (u  -V  ) 

'2 

^  =  ~~p  p  p  “  2(n-l)[E  +  F  +  (n-2)  [H  +  |(n-3)j]) 

2v  {m  -v)  ^ 

T  2^ 

h  P 

®  =  g  '3  P  ■  [D  +  F  +  “(n-2)(n-3)l] 

2uiu^-v^)  5 

A  =  -  2n(B  +  C  +  (n-l)(D  +  E  +  2F  +  i(n-2)[2H  +  (n-3)(l+J) ] }) 

For  the  n-cube,  w(X)  =  1;  u,v  =  [5/9[l  ±  (8/35)^^^])^^^  .  For  the  n-sphere, 

w(X)  =  l:  u,v  =  C5/(n+8)[l  ±[(2ri+6)/(5n+30)]^/^^/^  .  For  the  n-cube, 

x(X)  =  ]f  (l//l  -  (x^)^):  u.v  =  (5/8[1  ±  (1/5)^/®])^/®  .  For  infinite  n-space, 

n 

w(X)  s  exp(-^  (x^)^):  u,v  =  [5/2[l  ±  (2/5)^/^])^/^  . 
i=l 


Notice  that  we  can  obtain  two  precision  nine  integration  formulas 
by  Interchanging  the  role  of  u  and  v  in  the  solution  to  the  above  non¬ 
linear  algebraic  equations.-^ 


^  Actually  we  could  have  obtained  an  infinite  number  of  solutions  by 

putting  G  ■  kH,  k  an  arbitrary  constant.  In  this  case  we  could 
not  take  advantage  of  the  reduction  of  points  occurlng  as  a  result 
of  putting  G  ■  0. 
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The  numerical  values  of  the  coefficients  and  zeros  of  the  precision 
three  and  five  integration  formulas  are  readily  evaluated  by  computer.  In 
Table  V  (Appendix)  we  list  zeros  and  weights  of  particular  precision  nine 
integration  formulas  we  previously  discussed.  In  Appendix  A  we  give  examples 
to  illustrate  the  accuracy  of  integration  formulas  developed  in  this  thesis. 


Table  IV  is  an  excerpt  from  Stroud  [21]  to  which  we  have  added  our 
own  results.  This  Table  shows  the  results  of  some  approximate  calculations 
for  the  two  integrals 


(I) 

(II) 


x^x^x^^  j  Ij  2,  3,  4 

e  dx  dx  dx'^dx 


2  1 
^  -  2 


|x^- I  +  |x^-“  j)  dx^dx^dx^dx^ 


These  integrals  were  calculated  numerically  by  the  following  methods: 

METHOD  A.  The  integral  is  approximated  by 
m 

™  I 

j  =  l 

where**  =  {[yj2/2],  [^]»  [>/Io])  . 

The  calculations  for  Integral  I  were  given  by  Davis  and  Rabinowitz  (see  page 
29).  The  calculations  for  II  were  given  by  Stroud  [21]. 

METHOD  B.  The  formula  of  precision  2  given  by  Stroud  (see  page  I9).* 

METHOD  C.  The  formula  of  precision  3  given  by  Stroud  (see  page  lO). 

METHOD  D,  The  formula  of  precision  5  given  by  Htimmer  and  Stroud  [7]. 

METHOD  E.  The  Cartesian  product  of  the  two  point  formula  for  the  line  segment. 

METHOD  F  and  G,  The  formulas  of  precision  9  derived  in  this  thesis. 


**  The  symbol  [Y]  for  example  indicates  the  fractional  part  of  the 

number  Y. 

*  These  formulas  are  also  given  in  this  thesis. 
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TABLE  IV 

COMPARISON  OF  A  QUASI-MONTE  CARLO  METHOD  WITH  FIVE  QUADRATURE  FORMULAS 


METHOD 

m 

I  VALUES  OF 

INTEGRALS  II 

Exact 

;  value 

1.0693976 

1 . 0000000 

A 

4 

1.0646192 

.8523646 

8 

1.0592766 

. 9654735 

16 

1.0615567 

•9597948 

32 

1.0626119 

.9837211 

64 

1.0586261 

.9915821 

128 

1.0657314 

.9936179 

256 

1.0673119 

.9977914 

512 

1.0668403 

. 9998982 

1024 

1.0681499 

2048 

1.0685418 

4096 

1.0685545 

8I92 

1.0688021 

16384 

1.0691568 

32768 

1.0691964 

B 

5 

1.0686301 

1.0310313 

C 

8 

/ 

1.0622464 

.9855997 

D 

33 

1.0688327 

.8606629 

E 

16 

1.0693883 

1.1547005 

F 

177 

1.0694046 

. 9448504 

G 

177 

1.0693986 

. 9448505 

The  integrand  in  II  has  a  discontinuous  derivative  in  the  range  of 
Integration  and  is  therefore  weighted  heavily  against  the  formulas  derived 
in  this  thesis  and  the  Gaussian  formulas.  Hence,  although  the  8-point  form¬ 
ula  C  gives  better  results  than  the  33-point  formula  and  the  177-point 
formulas  this  will  not  generally  be  so. 

3)  Orthogonal  Polynomials 

By  our  method  of  setting  up  non-linear  algebraic  equations  we  obtain 
numerical  integration  formulas  of  precision  p  over  fully  symmetric  regions 
In  where  n  >  (p“l)/2.  We  shall  prove  the  validity  of  the  solution  to 
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the  set  of  equations  (9)  by  presenting  another  approach  to  the  problem  of 
obtaining  numerical  integration  formulas  over  fully  symmetric  regions  in  e'^; 
this  approach  will  prove  to  be  an  aid  in  finding  numerical  integration  form¬ 
ulas  of  precision  4k  +  1,  k  an  integer. 

i)  Polynomials  orthogonal  over  the  n-cube;  w  Ll) . 

As  an  aid  in  finding  integration  formulas  that  have  some  maximum 
precision,  we  construct  a  set  of  polynomials  orthogonal  over  the  n-cube 
according  to  the  following  rules*. 

(a)  All  the  Independent  variables  are  considered  equally  important.  For 

2 

example,  in  the  case  n  =  2,  if  the  orthogonal  polynomial  has  the  term  kx  , 
it  shall  also  have  the  term  ky  (i.e.  the  polynomials  are  symmetric  in  the 
variables ) . 

(b)  The  integral  of  the  product  of  a  monomial  degree  less  than  m  and  an 

orthogonal  polynomial  of  degree  m  over  the  region  shall  be  zero. 

By  these  rules  we  are  merely  trying  to  extend  the  idea  of  orthogonal  polynom¬ 
ials  in  one  variable. 

In  trying  to  construct  a  set  of  orthogonal  polynomials  in  higher 

dimensions  according  to  these  two  rules  we  find  that  there  exists  no  unique 

2  2 

set  of  polynomials.  For  example  each  of  the  bases  1,  x,  y,  x  ,  y  ,  ..., 

Ic  "Ic  lie 

x^,  y  ,  or  1,  xy,  ....  (xy)  ,  ...  could  be  used  to  construct  a  set  of 

polynomials  orthogonal  over  the  square.  Once  we  have  constructed  a  particu¬ 
lar  set  of  orthogonal  polynomials  we  may  moreover  find  that  the  zeros  of  each 
member  of  the  set  are  not  uniquely  determined,  since  they  are  level  curves. 
The  polynomial  x  +  y,  for  example  has  an  infinite  number  of  zeros,  x  =  -y. 

We  will  call  this  set  of  polynomials  g”(K);  C  for  cube,  m  being 
the  degree  and  n  the  dimension.  To  preserve  total  degree  the  pair 
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ii 

c'lix)  =  1,  c^(x)  =  V  is  an  obvious  one  to  start  with.  From  these  two 

i=l 

it  is  easy  to  see  that  all  other  orthogonal  members  of  the  set  will  have  mon¬ 
omials  with  either  odd  degree  or  even  degree,  but  not  both. 


n  i 


It  is  not  immediately  obvious  which  one  of  O^Cx)  =  1  +  a 


i  j 

X  x*^ , 


n 


=  1  +  (x^)^  ^  y  ’ 


or  1  +  a 


7 


i=l  j=l 


we  should  take.  We 


i=l 


i=l 


cannot  solve  for  b  in  the  second  case.  Consideration  of  C^(X)  tells  us 
nothing  we  don't  already  know,  so  let  us  consider  n  =  2.  Here 

we  have  the  following  choices; 

(i)  C^(X)  =  l+a(x^+y^+xy)+b(x^+xj+x^y^+xy^+y^) 

(ii)  C^(X)  =  l+a(x^+y^)+cxy+b(x^+y^)+ex^y^+d(x^y+xy^) 


(iii)  C^(X)  =  l+a(x^+y^)+b(x\y^)+ex^y^ 


(iv)  C^(X)  = 


l+a( x^+y^ )+b ( x^+y^ ) 


All  the  constants  can  be  uniquely  determined  in  only  the  last  of 
these  choices.  Hence  if  we  further  add  that  the  polynomials  should  be  the 
simplest  possible  that  satisfy  (a)  and  (b)  above,  we  find  these  polynomials 
to  be 


n  m 

\“ 


(10)  < 


i=ij=i 
n  m 

i-1  j=l 


By  simplest  possible  we  mean  that  the  polynomials  have  the  fewest  possible 
monomials  required  to  satisfy  rules  (a)  and  (b)  above. 


/ 


'■  .3. 

tiz;  ry/ 


■  ■>  .4:^  ,v  t  j , 


.  1 


wo  a; 
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%  fO'} 


i  ^■ 


:■:  is,;.-  'ij  sptfciq  ::; ds :•* i 8 

jiite.-,o::3  hauXupa^^i'  aIalm«^fi:oa:t 


By  evaluating 


(11)  (^J  ^  C"(X)  dX  =  0 


-1 

for  1  <  k  <  m  it  is  easy  to  show  that  the  coefficients  a^, 
c^(x)  can  be  obtained  by  solving  the  matrix  equation 


m 


of 


(12) 


11  1  ' 

r 

3  5  '  *  *  2m+l 

1 

11  1 

1 

5  7  ;  ‘  ’  2nH-3 

^2 

1 

3 

11  -1 

• 

n 

i 

2k+l  2k+3  “  •  •  2m+2k-l 

^k 

2k- 1 

11  -1 

• 

• 

1’ 

2ro+l  2m+3  ’  *  *  4m- 1 

Si 

m 

2iti-1_ 

Similarly,  the  coefficients 
solving  the  matrix  equation 


a  of  G_  1  (X)  can  be  obtained  by 


m 


(13) 


1  1  _ 1_ 

1 

5  7  ‘  '  2m+3 

3 

1  1  1 

a 

1 

7  9  ...  2nH.5 

2 

5 

• 

=  -  1 

ft 

1  1  *1 

« 

2k+3  2k+5  ‘  ’  ‘  2m+2k+l 

• 

\ 

2k+l 

• 

1  1  ’  •  1 

• 

• 

1 

2m+3  2m+5  .  .  •  4m+l 

a 

m 

2m4‘l 

Considering  the  case  n  =  1  in  equations  (12)  and  (I3)  the  poly¬ 
nomials  c”(X)  are  readily  Identified  with  the  Le^gendre  Polynomials  by 
the  equation 


n 

(14)  o“(X)  .  Y,  P„(xh 

l&l 


m 


where 


is  a  constant  of  proportionality. 
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It  is,  moreover,  easy  to  see  that  for  m  >  0, 


whenever  each  integer  k,  satisfies  0  <  k.  <  ra  . 

i  =  i 

ii)  The  choice  of  zeros  of  the  polynomials 

For  m  >  0,  n  >  1,  the  equation  =  0  has  an  infinite  number 

of  roots.  Moreover,  constructing  interpolation  polynomials  using 
presents  extreme  difficulties.  If,  however,  we  chose  the  zeros  in  a  manner 
corresponding  to  the  way  we  have  set  up  the  non-linear  algebraic  equations; 
that  is  by  first  setting  all  variables  except  one  equal  to  zero  and  then 
equating  the  resulting  polynomial  in  one  variable  to  zero,  we  find  these  zeros 
to  be  the  required  non-linear  unknowns  in  our  simultaneous  non-linear  alge¬ 
braic  equations.  For  example  the  zeros  of  C^(X)  chosen  in  this  way  are 
u  =  +  those  of  C^(X)  are  u  =  0,  +  (3/5)^^^,  and  those  of 

Cp(X)  are  u,v  =  0,  +  these  being  the  solutions  to  the 

non-linear  algebraic  equations  (7)>  (8)  and  (9)  respectively  for  the  n-cube„ 

The  zeros  of  chosen  as  above  will  always  be  the  zeros 

of  P_  i(X),  The  zeros  of  C^(X)  are  unfortunately  complex  for  m  >  1  -- 
it  can,  moreover,  be  shown  that  the  complex  zeros  of  C^(X)  are  not  the 
solution  to  the  non-linear  algebraic  equations  set  up  in  the  manner  described 
above  for  obtaining  an  integration  formula  of  precision  7- 

Hence  only  the  polynomials  useful  for  obtaining 

nauterical  integration  formulas  of  precision  p  i|k  +  1  for  n  >  (p*l)/2  . 

n 

lit)  Polynomials  orthogonal,  over  an  arbitrary  fttlly  symmetric  region  in  E 
Since  our  even  polynomials  ware  not  useful  for  the  specific  case. 


the  n-cube,  a  set  of  even  polynomials  will  not  be  useful  in  the  general  case. 


'  1  f'  1 


r,  ,■  I 
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The  odd  ones,  however,  prove  to  be  helpful.  We  define  these  to  be 

n  m 


i=l 


j=l 


i=l 


where  the  obtained  by  solving  a  system  of  m  simultaneous 

linear  algebraic  equations,  the  k'th  of  which  may  be  written 
m 

(17)  y  a  r  =  -  f  w(X)(x^)2^  dX  , 

L_J  J  Jr  Jr 

j=l  n  n 

and  w(X)  being  fully  symmetric  as  described  above. 


Although  these  polynomials  are  orthogonal  over  R  to  any  monomial 

.  k.  ^ 

i  \  1 

of  the  form  (x  )  ,  they  are  not,  in  general,  orthogonal  to  any  monomial  of 


«  n 

n  c  k , 

the  form  ||  (x^)  ^  where  0  <  /  ^*5 

Ul 


2m 


The  zeros  of  Q^,,(X)  chosen  in  the  manner  described  above  are 

i. 

obviously  those  of  i{x).  Considering  the  properties  of  R  and  w(i 

2m4-i  n 

we  may  write*,  with  the  x^  ®s  being  the  zeros  of  *^2m+l^^^’ 

m^ 


k 


r  „(X)  n 

1=1  J 


dX 


n 


j=l 


f  W(x)  n  (x"-x^^) 


dX 


Jr 


j=l 


-/n^  /2  n\ 

0_(x  )  p0„(x  ,...,x  ) 


-eaCn")  \(x2 . x")  ^=1  ^ 


-a 


Due  to  the  symmetry  of  the  region  R^,  the  functions  SgCx 
Q  (x^,...,x’^)  are  all  even  functions  of  their  independent  variables.  If,  in 


n 


addition  we  assume  these  functions  positive  throughout  R  (but  not  necessar- 

n 

ily  on  the  boundary)  we  may  write 


We  assume  here  that  R  contains  the  origin,  and  that  the  limits 

III 

of  integration  can  be  written  as  we  have  written  them. 
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r 


w 


wn 

j=i 


/  i  i  \  ^  ^  \ 

( X  -X  ,  (x  ) 


^a 


dX  =  j  ^(k  ) 


(lx 


-a 


2 

where  u)(x  )  is  positive  throughout  the  interval  (-a,a).  Hence  by  well 
known  arguments  (see,  for  example  [Ij])  all  the  zeros  of  ‘^2m+l^^^  real, 

distinct,  and  lie  in  the  open  interval  -a  <  x  <  a  . 


Again,  by  well-known  arguments,  the  zeros  of  the  polynomial  q  i(x) 
.  2rn-f  1 

be. 

can  be  shown  toy^the  non-linear  unknowns  in  the  non-linear  algebraic  equations 
set  up  in  the  manner  we  have  set  them  up  (e.g.  for  the  cases  p  =  5>  P  =  9) 
for  obtaining  integration  formulas  of  precision  p  =  4m  +  1,  n  >  (p-l)/2. 

Our  reasons  for  this  statement  are  as  follows:  Assuming  2m+l  distinct  non¬ 

linear  unknowns  (including  zero)  we  can  always  obtain  from  these  combinations 
of  fully  symmetric  sets  of  points  which  will  enable  us  to  integrate  all 
monomials  up  to  degree  4m+l  (although  we  have  not  yet  found  any  examples 
for  which  equation  (lO)  does  not  hold,  restrictive  conditions  such  as  these 
may  prevent  us  from  obtaining  solutions  in  all  cases).  The  2m  equations 
resulting  from  integrating  (x^)^,  (x^)^,  ...,  (x^)^  uniquely  determine  all 
our  non-linear  unknowns  (for  proof  see  [14]).  We  can  then  solve  the  result¬ 
ing  system  of  linear  algebraic  equations. 


2 

For  example,  if  we  multiply  equation  (ii)  of  (9)  by  u  ,  subtract 

2 

(iii)  from  it,  call  the  result  a,  (iii)  of  (9)  by  u  ,  subtract  (iv)  from 

2 

it,  call  the  result  b,  (iv)  of  (9)  by  u  ,  subtract  (v)  from  it,  call  the 
result  c,  then  v'^  =  a/b  =  h/c  is  a  polynomial  equation  whose  solution  is 
(11),  these  roots  being  the  same  as  those  of  q^(x)  (excluding  the  zero 
root).  This,  then,  is  a  partial  proof  of  the  validity  of  the  solution  to  the 
algebraic  equations  (9).  We  may  easily  obtain  the.  value  of  F  by  subtract¬ 
ing  equation  (viii)  from  (ix).  By  successively  eliminating  the  unknowns  I, 
J,  ...  from  equations  (vi),  (vii),  (viii),  (x),  (xi),  (xii)  we  can  solve 
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for  D,  E,  G,  H,  I,  J;  a  check  being  that  ~ 

two  of  equations  (ii)  and  (v)  will  give  the  solution  to  B  and  C.  From 

equation  (i)  we  can  then  obtain  the  value  of  A. 

In  carrying  out  this  outlined  procedure  we  would  prove  the  validity 
of  the  given  solution  to  the  non-linear  algebraic  equations  (9)* 

iv)  Polynomials  orthogonal  over  the  n-sphere 

For  completeness  we  also  find  the  polynomials  ~  ^^+1^^^ 

orthogonal  over  the  n-sphere  with  respect  to  the  weight  function  1. 


We  transform  the  integral  over  the  n-sphere  to  an  integral  over  an 
n  dimensional  parallelepipedon  using  polar  coordinates.  Now  for  the  circle 
(2-sphere)  the  transformation  that  will  do  this  for  us  is 


(18) 


1  „1 

X  =  rcos0 

2  .  ^1 
X  =  rsinO 


(19) 


Similarly  for  the  sphere  (j-sphere)  we  may  let 


r  1 

X  =  rcos0 

2  12 
X  =  rsin0  cos0 

3  ,1,2 

X  =  rsin0  sin0 


to  transform  an  integral  over  a  sphere  to  an  integral  over  a  rectangular 
region  in  three  dimension.  Both  (18)  and  (I9)  satisfy 


n 


(20) 


/  iv2  2 
(x  )  =  r 


i=l 


and  hence  using  (20)  in  view  of  (18)  and  (I9)  the  transformation 


-  56  - 


(21)  < 


1  „1 

X  =  rcosO 

2  12 

X  =  rsin0  cos0 

5  12^ 

x"^  =  rsinO  sin0  cosO"^ 


^  =  rsin0^sin0^.  .  . sin0”  ^cos0’^  ^ 

n-1  .-1.-2  ,  -n-3  .  ^n-2  ^n-1 

X  =  rsin0  sin0  ...sin0  '^sln0  cos0 

n  ,-1.-2  .  -n-3  .  ^n-2  .  .,n-l 

X  =  rsin0  sin0  ...sin0  '^sin0  sin0 


extends  (l8)  and  (I9)  to  the  n-sphere,  where  n  >  2. 

We  next  proceed  to  find  the  Jacobian  of  the  transformation 

(21).  Here  it  becomes  convenient  to  define  the  more  concise  notations 


(22) 


^  k  k 

p  =  COS0 

k  .  k 
V  =  sin0 


J  =  f 
n  /  n' 


n-1 


Then  from  the  well-known  definition  of  the  Jacobian  =  | j j  where  on  trans¬ 


forming  from  the  variables  x  to  the  variables  x,  a 


3X 


IJ 


it  is 


easily  shown  that  in  terms  of  our  definition  (22)  and  equation  (21) 


(23)  J  =  (see  page  57). 

n 


On  inspecting  the  last  two  rows  of  this  determinant  we  note  that 

1=1  1=1 
//  n-lv2  /  n-lv2v.,  ^fr  i 

=  ((1^  )  +  (v  )  )J^_^  11  V  » 


(24) 


i=l 


n-2 


J  =  J  H 
n  n-1  iJ- 


i=l 
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•H 

> 
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1 
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2 

2 
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since  =  1.  Since  J  =  1,  J  =  v^,  we  have 

lNn-2/  2vn-3  ,  n-3v2/  n-2 


J  =  (v‘)“  "(v") 


(25)  /;  =  r"-^  n 

^  1=1 


(v  (v  ),  from  which  we  deduce  that 
n-2 

n-i-lvi 


For  purposes  of  integration  over  the  n-dimensional  unit  sphere  the 
ranges  of  variables  may  be  chosen  to  be 

0  <  r  <  1 


(26)  0  <  0  <  «  1=1,... ,n-2 

0  <  <  2it 


and  we  have  now  arrived  at  the  result  that 


(27) 


r...  r  f(x^, . , , ,x”)dx^. . .dx° 

J  \.J 


^  I  *[{  ((sin0""^"^)^d0'-)d0 

0  0  0 

where  F(r ,0^, . . .  j©*'  =  f  (rcos0^, . . .  ,rsin0^sln0^. . . sln©’^  under  the 

transformation  (21). 


n-2 


n-2 


(28) 


In  terms  of  the  transformation  (21)  and  the  notations  (22)  we  have 
n 

R  W  r  k  It  k  k 

n  /  i\  i  1/  In  1/  1  2x  2  /  1 ,2  n-2  n-lx  n-1 

[  (x  )  =  r  (|i  )  (v  M  )  ...(v  V  ...V  ) 

!■  i  u 

t  1  2  n-l\  n 
•(v  V  ...V  ) 


n-1  ,  k,  ,  ?  ,k 

^  n  ((t'  )  hvh  ^ 

i»l 


Substituting  k  for  n-i-1  in  (25)  and  writing  1  for  k  yields 

(£9)  j„  ■  "n  .  ”n  . 

"  1-1  m 


4' 


'I  /  \  ■  /  ;  ^  4. 


.1 


‘•‘V  ^ 


<..}  x  sn  j  1 


.  »  _  ^ 

’ll  1.  ^ : 


h 
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On  multiplying  the  right  hnnd  side  of  (28)  by  the  right  hand  side 
i  i  i  i 

of  (29)  and  writing  cos0  for  ,  sin0  for  v  we  obtain  the  integ¬ 
ration  formula 

(30)  f ■■■  f  n  (='9  ^ 

^  J  1=1 

n-sphere 


n. 

L-!-  E  k 
s=i+l 


0  0 


If  at  least  one  of  the  k^*s  is  an  odd  integer  the  integral  (50)  is  zero  due 
to  the  symmetry  of  the  region  of  integration.  If,  however  all  the  k^*s  are 
even  positive  integers  or  zero,  (30)  will  be  different  from  zero.  To  evalu¬ 
ate  (30)  we  borrow  a  well-known  result  from  the  theory  of  the  Gamma  function, 
namely 


(31) 


(sine)kpose)-'de  =  |  ^ 


H..  1  r(!^)  r(J|l) 


where  in  our  applications  we  assume  j  and  k  are  zero  or  positive  integers 


Writing  2k for  k^  in  (30)  and  using  (3I)  we  obtain 


^  2k 


/•••/  lOi 


^  dx^ 


n-sphere 


n 


g  "-i  r(-|-)  r(-— 


.E,(2k.+1)  2k.4-n-i+_  ?^i2k 

iffll'  i  '  1=.  ^  1  _  sai+i 


Due  to  cancellation  of  successive  terms  in  the  above  product,  (32)  reduces  to 

n 

1  .  2k, 
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n-sphere  * 


“i  j  1 
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which  is  in  agreement  (by  putting  each  =  O)  with  the  well-known  result 
for  the  hypervolume  of  the  n-sphere* 


n/2 


(5^)  V  =  r. . .  r  n  dx^  =  — 

n  j  j  -- « 

n-sphere 


r(i+|) 


Using  formula  (5^)  together  with  the  notation  =  a(a+l) . .  .  (a+k-l) , 

(a)^  =  1,  and  the  recurrence  relation  r(l+x)  =  xr(x)  for  the  Gamma  func¬ 
tion  we  finally  obtain 


(55) 


n 

r.r  ii 

1  ^2  ''(.Z,  k, ) 

n-sphere  'i^rl  i' 


.  V 


n 


We  can  now  write  down  the  set  of  linear  algebraic  equations  from 

which  we  can  find  the  values  of  the  coefficients  a,  .....a  of  n(X)  = 

1  m  2m+l'  ' 

t(X).  For  on  substituting  the  appropriate  form  of  (35)  into  (I7)  and 
simplifying,  we  find  these  equations  to  be 
(2r-<-lv 

(36)  )  ^  2^3  .a  =  -1,  r  =  1  to  m. 

Equations  (35)  and  (36)  can  be  used  to; 

(i)  find  the  explicit  solution  to  equations  (9) 

(il)  verify  the  solution  to  equations  (8) 

(iil)  verify  that  (lo)  holds 

(iv)  find  numerical  integration  formulas  of  higher  precisions  for 
the  n-sphere. 

n 

If  we  write  S^^^(X)  =  ^  then  according  to  the 

Inl 

result  of  pages  53“54  all  the  aeros  of  Sg^^i(x)  lie  in  the  open  interval 

-1  <' X  <  1, 
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CHAPTER  IV 

INTEGRATION  FORMULAS  OF  ARBITRARY  PRECISION 


In  Chapter  III  we  have  developed  a  relatively  simple  method  of 
obtaining  numerical  integration  formulas  by  reducing  the  solution  of  a  system 
of  non-linear  algebraic  equations  to  the  solution  of  a  system  of  linear  alge¬ 
braic  equations.  The  formulas  we  can  obtain  by  this  procedure  are  of  prec¬ 
ision  p  over  fully  symmetric  regions  in  n-space  where  p  =  Uk  +  1,  k  an 
integer.  To  overcome  the  restriction  this  equation  imposes  on  p  we  call 
on  the  repeated  Gaussian  formulas,  noting  (see  Table  II  for  example)  that  for 
large  p  and  moderate  n  these  formulas  are  not  unduly  extravagant  in  ordin¬ 
ate  evaluation,  particularly  since  they  yield  greater  accuracy  for  larger  values 
of  p . 


Arbitrary  limits  of  integration  in  an  n-fold  integral  are  somewhat 
inconvenient  and  the  practical  task  of  evaluation  is  eased  if  we  can  transform 
the  limits  to  (-1,1),  i.e.  if  we  can  transform  an  integral  into  an  integral 
over  the  n-cube.  A  transformation  is  given  which  enables  us  to  do  this. 

In  this  chapter,  formulas  of  arbitrarily  high  precision  over  the 
finite  and  infinite  n-sphere  are  constructed.  These  formulas  are  included 
for  completeness  rather  than  for  practical  value.  In  practice,  we  are  re¬ 
stricted  to  values  of  p  and  n  which  can  be  encompassed  in  a  reasonable 
time  on  existing  hardware.  The  formulas  can  however  achieve  remarkable 
accuracy  in  evaluating  integrals. 


l)  Formulas  over  Rectangular  Regions 

The  methods  of  Chapter  III  suffice  for  formulas  of  precision 
p  =  1+k  +  1,  k  an  integer. 


It  is  desirable  that  formulas  of  all  odd  precisions  be  available, 
particularly  for  moderate  values  of  n.  The  main  difficulty  impeding  the 
development  of  economical  high  precision  formulas  is  the  brute  difficulty 
of  solving  a  large  system  of  non-linear  equations  in  a  number  of  unknowns. 
No  general  methods  are  available  and  a  search  of  the  literature  suggests 
that  little  research  is  being  done  in  non-linear  systems.  In  this  Chapter 
we  evade  the  difficulty  by  using  Gauss-product  formulas. 

Suppose  the  region  of  integration  is  a  rectangular  region  in 
n-space.  The  Gaussian  integration  formulas  produce  very  good  accuracy  for 
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the  large  class  of  functions  that  have  all  their  derivatives  continuous  in 
this  region  of  integration.  Moreover  the  number  of  evaluation  points  re¬ 
quired  is  for  fixed  n  the  number  of  points  increases  polynom- 

ially  in  p,  not  exponentially.  Thus  if  the  number  of  dimensions  in  n- 
space  is  small  (e.g.  2,  >  5)  we  may  be  able  to  obtain  accurate  results 

with  a  feasible  number  of  points. 

The  construction  of  product-type  integration  formulas  over  rect¬ 
angular  regions  in  n-space  is  not  difficult.  Let  the  region  of  integration 
be  described  by 

i  i  ,  i 
a  <  X  <  b  -x-* 

and  let  the  weight  function  w(X)  have  the  form 

n 

w(X)  =  ]]  w  (x^)  , 

i=l 

1  i  i  i 

w^(x  )  >  0  in  a  <  X  <  b  ,  then  if  we  can  find  polynomials  orthogonal 

over  the  intervals  (a^,  b^)  with  respect  to  the  corresponding  weight 

functions  Wj,(x^)  we  can  find  the  Gaussian  integration  formulas  of  arbit- 

i  i 

rary  precision  over  the  intervals  (a  ,  b  )  with  respeset  to  the  weight 
functions  w^(x^).  Using  the  Cartesian  product  theorem  given  on  page  11 
we  easily  obtain  numerical  integration  formulas  of  arbitrary  precision  over 


* 


•x* 


We  actually  have  a  variation  here;  more  generally  we  require 
TT  Pi'^^ 

li  ('"2""")  where  p^  is  the  precision  obtained  in  each  variable. 


1  i 

Here  a'  and  b  may  be  finite  or  infinite,  assuming  the 

weight  functions  w^(x  )  to  be  appropriately  chosen. 
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rectangular  regions  in  n-space. 


n 

If  it  is  not  possible  to  write  w(X)  =  |T  w  (x^),  or  if  some 

i=l  ^ 

of  the  w^(x^)  ^  0  everywhere  in  (a^,  b^),  we  can  always  include  such 
weight  functions  as  part  of  our  integrand  and  apply  standard  Gaussian 
formulas  (Legendre,  Laguerre  and  Hermite);  e.g.  if  we  define  w(X)  f(X)  = 
F(X),  we  have 


(1) 


a 


pb 

/  w(X)  f(X)  dX  = 
J  n 

a 


For  example  if  the  limits  of  integration  are  finite  we  know  that  we  can 
employ  the  Gauss -Legendre  formulas  to  evaluate  the  above  integral  on  the 
right  to  arbitrary  precision*.  Due  to  the  additional  non-linearity  of  the 
weight  function  w(X)  we  will  not  in  general  be  integrating  f(X)  to 
precisions  p^  (see  footnote  on  previous  page)  when  integrating  F(X)  to 
precisions  p^  , 


2)  Transformation  from  an  Arbitrary  Region  to  the  n-Cube 

The  results  of  this  section  make  it  possible  to  apply  the  repeated 
Gaussian  integration  formulas  to  integrals  of  the  type  (2)  below.  The 
following  theorem  enables  us  to  carry  out  the  transformation. 

THEOREM.  Given  the  integral 


* 


If  some  of  the  limits  of  integration  are  infinite  we  can  employ 
Gauss -Laguerre  integration  (and  similarly,  Gauss-Hermite)  by 
setting  P( ■ )  "  exp(-x*)  0(*)!  0(*)  «  exp(x*)  F(*)  . 
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(2) 


,1  ,  2/  1  X  ;  1  11-1  V 

(y.)  (x . X  )  2  n,,  n,  b-1  ,  1 

I  =  /  /  ...  f  f  (x  ,x  ,  .  .  .  ,x  )dx  die  ...dx 

U  Jn.  1  n-lx 

qi  cp(x)  cp(x,..,,x  ) 


we  can  always  transform  this  into  an  integral  over  the  n-cube  by  a  sequence 

of  n  linear  transformations. 

i  i  11 

PROOF.  The  symmetric  case  cp  =  ,  i  =  1 , .  .  .  ,n  (cp  and  ijj  are  constants) 

is  easily  tractable  and  we  consider  it  first.  We  set 


(3) 


i  ,i  i 
X  =  ij/  u 


i  =  1, . . . ,n  . 


Then 


dx^ 

dx^ 

ih) 

• 

= 

J  ^ 

dx 

where  by 

,  i 

we  mean 

r 

2,2 

u 


0 

0 


n.n  n.n  ,n 

uip^  u  ^2  •  •  •  ^ 

.i 


du‘ 


du 


du 


n 


,  From  (4),  the  jacobian  J  of  the  transformation 
(5)  is  j  ust  the  product  of  the  diagonal  elements  of  the  matrix; 


(5) 


J  = 


i 


Hence  we  have 


(6) 


,1  iH,  1  n-lv 

ri"  . *  ^  1  n,  .  n  ,  1 

I  •  •  •  /  f (x  , . , . ,x  )  dx  .  .  .  dx 

1  .n,  1  n-lx 

(x  , .  .  .  ,x  ) 

r  ...  r  F(u^, .  . .  ,u’^)  J  du’^...du^ 

d  J 

-1  -1 


where  J  is  given  by  (5),  and  F(u^,...,u”)  =  f  (x^ , . . .  jx’^)  under  the  trans¬ 
formation  (5). 


If  the  limits  of  Integration  are  as  in  (2)  then  the  set  of  trans¬ 


formations 
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(7) 


i  iIj 

X  =  ^ 


±JL 


+  u 


2 


i  =  1, . . .  ,n 


will  transform  the  Integral  (2)  to  an  integral  over  the  n-cube,  and  the  proof 
is  similar  to  that  for  the  symmetric  case  above.  The  Jacobian  of  the  trans¬ 
formation  (7)  is 

(8)  j  =  . 


This  theorem  thus  allows  us  to  apply  Gaussian  integration  formulas 
to  integrals  of  the  type  (2).  It  may  be  worth  while  finding  integration 
formulas  with  the  functions  J  as  weight  functions  if  J  is  suitable  and 
the  value  of  integrals  are  required  over  the  region  for  a  sufficient  number 
of  functions  f(’)*  require  the  value  of  only  one  integral  it  will 

likely  be  more  economical  of  programming  time  to  use  standard  Gaussian  inte¬ 
gration  formulas  (e.g.  Legendre,  Chebychev)  applicable  over  the  interval 
(-1,1)  with  more  points  to  sustain  accuracy.  Considerable  extra  expense  of 
machine  time  is  tolerable  for  one-shot  evaluates. 


Let  us  consider  the  simple  case  in  which  the  function  f(')  is  a 
polynomial  in  n  variables  x^(i=l, . . . ,n) .  Then  F( “ )  =  f(*)  under  the 
transformation  (7)  will  in  general  not  be  a  polynomial  in  the  variables 

1.  XI 

u  (i  =  l,...,n-l)  although  it  will  be  a  polynomial  in  u  of  the  same  degree 

as  f(*)  was  in  x’^.  Hence  the  transformation  (7)  will  introduce  complica¬ 
tions  in  n-1  of  the  variables  of  the  integrand  and  when  performing  numeric¬ 
al  Integration  tl.e  precision  obtained  ruay  not  be  as  high  in  the  origiT'^al 
variables  as  in  the  new  variables. 


As  an  example,  let  us  transform  the  integral  over  the  circle 
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onto  the  correspond ini^;  jntO£;ral  over  the  square.  By  equation  (i|)  the  linear 
transformations  which  will  enable  us  to  do  this  are  x  =  u,  y  =  = 

Wl-u^.  By  equation  (5)  the  Jacobian  of  the  transformation  is  L “ u^ ,  Hence 


VJe  can  now  perform  repeated  Gaussian  integration  to  evaluate  this  integral; 
Gauss -Lengendre  in  the  variable  v,  Gauss ian-Chebychev*  in  the  variable  u. 
The  resulting  integrand  f(u,v7T-u^)  is  however,  more  complicated  in  u  than 
f(x,y)  was  in  x  (even  though  the  Gauss -Legendre  integration  will  remove 
all  odd  powers  of  s/l-u®  in  f (u, w/l-^^)  . 

This  transformation  theorem  is  of  course  limited  in  scope  to  inte¬ 
grals  of  the  type  employed  in  the  statement  of  the  theorem.  Within  this 
limitation,  the  theorem  is  of  practical  value  since  it  enables  us  to  pres¬ 
cind  from  the  task  of  finding  integration  formulas  valid  for  arbitrary 
(bounded)  regions  in  n-space  and  to  concentrate  on  the  simpler  problem  of 
finding  higher  precision  integration  formulas  for  the  n-cube,  keeping  the 
weight  function  arbitrary. 

3)  Formulas  of  Arbitrary  Precision  over  the  n-Sphere 

This  problem  has  already  been  solved  by  W.  H.  Peirce  for  the  cases 
n  =  2  and  n  =  5  (see  [I5]  and  [16]).  The  formulas  developed  (here)  in 


* 


We  can  employ  the  integration  formula 

1  ^ 
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cos 
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The  X  's  are  zeros  of  the  Chebychev  polynomial  of  the  second  kind; 

J  . 

U  (x),  orthogonal  over  the  Interval  (-1*1)  with  respect  to  the 

weight  function  vl-x^  . 
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section  4  of  the  last  Chapter  enable  us  to  extend  the  results  of  W.  H.  Peirce 
to  the  n- sphere,  where  n  >  2. 


We  shall  seek  numerical  integration  formulas  of  a  fairly  general 
form,  namely 


(9)  S  S  f(x^,  .  .  .  ,x”)dx^dx^.  .  .dx’^ 

n-sphere 


m^  m,  m  , 

0  1  n-1 


-A 


£/  1  n  V  . 

a,  .  ,  fix,  .  .  , . . .  ,x,  ,  ,  )  , 

Vi-'-Vi  Vr-'Vi 


the  radius  of  the  sphere  may  be  finite  or  infinite  depending  of  the  weight 
function  w( * ) . 


(10) 


By  formulas  (27)  and  (29)  of  Chapter  III  we  may  write 

J' ' ' '  J'  )  f  (x^, . . .  ,x’^)dx^. . .  dx’^ 


n-sphere 


=  f  w(r)r"-l  n  (/”{sinel)"-l-l)  f 

0 


2jt 


F(r,0\0^,...,0’^"b 


. d0""^. . .d0^dr 


where 


F(r  ,0^, . .  .  ,0’^  =  f  (rcos0^, . .  .  ,rsin0^. .  . sin0’^ 


under  the  transformation  (21)  of  Chapter  III,  For  n  =  2  we  exclude  the 


product  ||(  P •)  in  formula  (lo). 


J 

Using  the  Cartesian  product  theorem  (page  ll)  and  the  transfor¬ 
mation  theorem  (page  9)  we  know  that  we  can  find  the  integration  formulas 
(9)  if  we  find  the  integration  formulas 


-1 
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(12) 


(13) 


111. 


r 


j 
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(sinO) 


n-k-1 


f,  (0)d0  = 

K 


/ 


£,  (0.  : 

Jk  Jk 


1<  =  1  Lo 


R 

r  w(r)  r”  ^  f  (r)dr  ^  V  c°  f  (r,  ) 

J  o  A  j  o  J 

0  i  =1 


Moreover,  in  view  of  the  transformation  (21)  of  Chapter  III,  if  the  formulas 
(13),  (12)  and  (11)  will  be  of  precision  p  in  the  variables  r,  sin0  and 
COS0,  formula  (9)  will  in  turn  be  of  precision  p  in  the  variables  . 


(14) 


Putting  y  =  cos(0/2)  in  (ll)  we  obtain 

dy 


J 

0 


- 1  - 1 

p2jt  pi  f  i(2cos  y)  A 

/  f  {9)dB  =  2  I  - T — ^ — ~  dy  ^  ^  c.  f(2cos  y.  ) 

'  >1  (l-y2)4  ,  Vl  J„-l 


this  being  the  well-known  Chebychev  integration  formula  (formula  (8)  Chapter 
l).  Chosing  =  2m  we  obtain  precision  p  =  km-l  in  the  variable 

y  =  cos (0/2);  or  p  =  2m-l  in  the  variable  y^  --  in  this  case  we  would 
also  obtain  precision  p  =  2m-l  in  the  variables  cos0  and  sin0.  The 
zeros  y,  and  weights  c,  for  formula  (14)  are  given  in  Table  VI. 


(15) 


Letting  y  =  cos0  in  (12),  we  obtain 


P  (sin0)”"^'^  f^(0)d0  =  J  (l-y2)(^“^“2)/2  f^(cos"^y)dy 

'0  -1 


f  (cos 


) 
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For  n  =  2  these  formulas  are  missing. 
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Here  we  assume  n  >  2,  k  =  1  to  n-2.  If  we  take  m  =•.  m  in  (15)  "se 

K 

the  zeros  and  weights  as  indicated  in  Table  VI  we  obtain  precision  p  =  2m-l 
in  y  =  COS0,  and  on  inspecting  the  transformation  (21)  of  Chapter  III  we 
note  that  the  errors  in  the  odd  powers  of  sin0  in  formula  (15)  will  cancel 
out  in  our  final  integration  formula  (9)  due  to  all  integration  points  being 
spaced  symmetrically  about  0  =  jt/2. 

We  consider  three  cases  of  the  integration  formulas  (I3)* 

(i)  w(r)  =  1  ,  R  =  1  ; 

(ii)  w(r)  =  e”’^  ,  R  =  00  ; 

2 

(iii)  w(r)  =  e  ^  ,  R  =  00 

We  assume  that  n  >  2. 

In  the  cases  (i)  and  (ii)  the  zeros  and  weights  required  are 
readily  obtained  from  the  references  indicated  in  Table  VI,  With  ra  points 
in  formula  (I3)  we  obtain  polynomial  precision  2m-l  in  r.  Case  (iii)  is 
readily  identified  with  case  (ii)  by  setting  y  =  r®  . 

We  have  used  the  notation  of  [18]  in  Table  VI. 


Hence  we  hay®  found  the  integration  formulas  of  arbitrary  precision 
over  the  finite  and  infinite  n~ sphere.  Formula  (9)  now  be  written 
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where  the  constants  c.  ,  k  =  0  to  n  -  1  can  be  obtained  from  Table  VI„ 

The  points  x.  .  can  be  obtained  from  the  references  indicated  in 

^0*  "-^n-l 

Table  VI  together  with  the  formulas  of  the  transformation  (2l)  of  Chapter 
III. 


We  see  that  the  total  number  of  points  required  to  obtain  this 
precision  (2m-l)  is  2(m)^.  This  number  of  points  (for  moderate  p  and 
large  n)  is  much  too  large  since  by  the  method  of  the  previous  chapter, 
formulas  of  precision  p  require  only  ©[n^'^  ^)/^]  points.  However,  the 
formulas  of  this  Chapter  are  not  too  extravagant  in  points  for  moderate  n 
and  large  p;  in  particular  we  have  integration  formulas  of  precision 
1  --  these  we  could  not  obtain  by  the  methods  of  Chapter  III. 
Further,  all  evaluation  points  lie  within  the  n-sphere. 
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CHAPTER  V 

ERROR  BOUNDS  FOR  n-DIMENSIONAL  INTEGRATION  FORMULAS  OF  GAUSSIAN  'TYPE 


To  obtain  a  bound  on  tlie  error  of  a  repeated  Gaussian  integration, 
we  proceed  as  McNaraee  has  done  in  [I9]  in  the  one-dimensional  case--in  this 
paper  he  expresses  the  error  of  Gaussian  integration  as  a  complex  integral 
and  then  proceeds  to  bound  this  error  on  a  suitable  contour  in  the  complex 
plane.  A  bound  on  the  value  of  a  complex  integrand  is  often  more  easily 
obtained  than  a  bound  on  the  derivative  of  a  moderate  or  high  order  of  a 
real  integrand. 

In  the  main,  the  analysis  is  restricted  to  integrands  which  are 
functions  of  two  variables;  corresponding  results  for  n-variable  integrands 
have  been  derived  and  are  stated  in  this  Chapter  without  proof. 

In  section  l)  of  this  Chapter  we  consider  Gauss-Legendre  integra¬ 
tion,  In  section  2)  we  present  the  analysis  of  the  repeated  Gauss-Leguerre 
case  and  write  down  the  corresponding  result  for  the  Gauss-Hermite  case. 

The  results  are  then  extended  to  repeated  integrals  that  are  combinations  of 
Legendre,  Laguerre  and  Hermite  integrations. 


1)  Gauss-Legendre  Case 

Denoting  the  zeros  of  the  Legendre  polynomial  P^(x), 

X.,  y,  respectively,  and  the  corresponding  weights  by  a  ,  b.  respectively, 
J  K  J  R 

we  can  obtain  a  bound  on  the  difference 


(1) 


^2  =  f  I  dx  dy  -  y  y  a^b^  f(Xj»y|^) 

'-1  '-1 


n  m 


j=l  k=l 


by  employing  contour  integration  and  asymptotic  techniques, 


Suppose  that  f(u,v)  is  an  analytic  function  of  two  independent 
complex  variables,  regular  within  and  on  the  boundary  of  regions  that  enclose 
the  strips  -1  <  Ri  u  <  1,  -1  <  R/  v  <  1  of  their  respective  complex  planes. 

Let  C  be  a  simple  contour  that  encloses  the  strip  -1  <  Ri  u  <  1  and  lies 
within  its  region  of  regularity  of  the  complex  u-plane.  We  can  then  deduce 
from  the  contour  integral 


WAIc.aUAt)  HO  CAJUi'iSOH  flOXTAJ^^SlTKi  JAi?0I8M3MXG"*a  iiJOH  StJP^lTOS  HOH&'S 


noi4-^tr^33fix  nBisRirfiO  ba.icvega-j;  ^  jo  ’xo'S’s^  orfj  no  JjfifciocF  &  nifi:ido  oT 
ai.d^  iftifoxaft3/?iX/>-“90o  sj;^f.l  rti  nj;  snob  eurf  ©©uifiWoM  sjj  baaoo'jq  3w 

Iftigfatjni  jci-ugwoD  .fi  ais  noidiiTgoJ'ax  qtiiaawBO  5o  *xo?is»  arf3  asees'iqjtja  ftrf  trsqfiq 
x«>iQmq,a  al  xrodfroa  a  no  joi'ts  ©Irfd  bquod  od  atj^^saotq  rmrfa  bfrr, 

530m  jlrodlo  ar  brix-igs^ril  Xolgwoo  m  io  guXfiv  9ff3  <io  bnuod  A  ,3fr;:fq 
&  io  3sfc30  rfgiiff  ':ro  o^fi^sboar  b  io  avi3;;vi35b  adi  no  bnued  b  jrcria  boniajdo 

.bnS3g!S3lRtX  |b03 

33B  Xioxdw  eferfB3go3fil  04  b93oX34eo3  ax  ele^^Iana  sri3  &d:3  nl 

8bqB'x394n.t  sX.dBxtev^n  adi  .'jditiaax  ^aJ-baoqfas j3go  ;8aXdBi3Bv  owd  io  anoldonrxi 
.io03q  dnoifilw  t5.iqBdO  axdi  fU  bsinaa  97®  hn&  baviiab  nsad  avarf 

-Bagair'i  oxbn9gs.J~eai.niT)  Tafaienoo  ow  xasqadJ  »XW4  io  (X  noiJoae  al 
9T-X3r;j9j-e3nf,0  bs^Bsqax  arf;i  io  eleYi»^f6  dnasaxg  sw  (2  nolioss  nX  .nois 
.asBo  s^xmx^H-eenfiD  odd  xoi  dXi;««'x  gfiXbfloqeaixoo  sdi  nwob  siiiw  fens  asao 
o  dfiOXdB.o.tdf/fOo  93B  dfidd  eXBt?g9inx  hoifisqai  o4  bobnssxB  nadd  ©7s  adXueax  ©dx 

.snoX3(i3gsdnX  ©dljxixoH  bn*  ©TtsitgsJ  ^.oxbaagdJ 


©grO  ,sxbfl©g9j«gBiisO  ( X 
IalfflO0^ioq  9ibxi©3©J  »d4  ao  rsoias  add  gnidonaO 
[©''’■tioaqesx  ».©  \;d  edrigiaw  Sfii.bnog89a7oo  add  bwa  «t^avjt3o»qQ©7  ^  ,  x 

sor.5xaiiXfe  «rf4  no  bnwod  s  ctXaddo  nao  aw 
m  XX 

J  '\  'I  "gS  (I) 

i-t  1-  X- 

.aaupXndoaa  oidodgm^^ee  bfta  rto.X36tg©dxU:  iwodnoo  gnixoXqm©  ’<d 


dnsbnsqsbnX  owd  3o  noxdonui  oid^^XanB  n©  aX  (v<w)2;  dadd  ssoqqnS 

iiCtro  dsdJ  &.aoXg97  io  YXBbnwod  add  no  bo*  oXddXw  XBlwgax  .aaXdBldBV  xaXqwosi 

t&iq  x©.tqmoo  sv'-rdoaqao-x  xxodd  io  X  >  v  4^3  >  I-  ,X  >  y  U  >  X-  aqi^de  add 

X  bn*  X  >  n  5.XI  >  I-  qXxda  odi  essoXona  ;vodd  dwpJrros  a  ad  0  dsJ 

o* 

(onn^b  fii>dd  rtoo  oW  .OffBiq.n  XaXqmoo  odj  io  qdXdAXngnd  io  noXga?  «dl  oXddiw 

XB'^x^gdnX  dwodfioo  add  ©oii 


-  73  “ 


^  r 

f(u,y)  du 

2ni  J 

C 

u 

(u  -  it)P^Ju) 

that 

m 

(2) 

Hi 

X 

II 

I  { 

f(x,,y) 


(x  -  X 


j=l 


111:  J. 

)P' (x.)  2jt 

m  J 


r 


m 


2jti  J  (u  -  x)P  (u) 
C 


du 


If  we  integrate  this  equation  with  respect  to  x  over  (“-1,1)  and  inter¬ 
change  the  order  of  integration  in  the  repeated  integral  on  the  right,  we 
obtain 

f(x.,y)  pi  P„(x) 
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pi  f(x  y)  pi  P  (X 
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dx 


+  -h  ^  IW  dxdu 

2ni  J  ?Jm)  j  u  -  X 
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We  can  now  employ  two  known  results.  If  u  is  not  a  real  number 
between  (-1,1)  —  unless  it  be  a  zero  of  P  (x)  then 
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2Q, 


t  \  p  (x) 

(u)  5=  /  m^  ' 

\.J 


dx 


u  -  X 


where  1®  Legendre  function  of  the  second  kind.  We  also  have 

Q  (u)  a  ™  P  (u)  log(^^^"4-)  “  f  ,(u) 

'  2  m  '  "'u  -  1'  m-1'  ' 

where  ®  polynomial  of  order  m  -  I  (the  coefficients  of 

are  tabulated  in  Whittaker  and  Watson  [20]). 


Using  these  results  we  can  write  (3)  as 
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The  vvoights  a  ,  are  defined  by 

pi  r  (x)  dx  2f  i(x.) 

/  m  * _ _  m- 1  ^  1  ^ 

~  J  (x-x  )p' (x  )  “  P'  (x. ) 

J  j  in  J 

they  have  been  extensively  tabulated  [29,  ^0,  Jl] 


Setting 

I  =  r  P  f(x,y)  dx  dy 
J  J 

-1  -1 


m 


we  have 

I  =  -j^y  a^f(x^,y)  +  Ej(y)j-  dy 

^-1  KL 

where  E|(y)  is  given  on  the  right  side  of  (4).  Repeating  the  process, 
have  in  view  of  equation  (1|) 


v;e 


m  n 


I  -  r  E^(y)  -Jy  =  £  ^  ^  ^2 

'-1  >1  kTi 


where 


(5) 


"2  = 


1  r  i=i 

—  /  2 -  Q  (v) 

■Ki  J  p^(v) 


dv 


being  an  arbitrary  contour  within  the  region  of  regularity  of  f(xsv), 
that  encloses  the  strip  -1  <  Ri  v  <  1  of  the  complex  v-plane.  Using 
equation  (4)  again  we  can  replace  the  sum  in  (equation  (5))  by  its 

corresponding  real  and  complex  integral: 


P  f(x,v)dx- 
J  J 

-1  C 


I  p  f(u,v)  Q^(u) 


du 


Recall  that  is  an  arbitrary  simple  contour  enclosing  the 

strip  -1  <  Ri  V  <  1  in  the  region  of  regularity  of  f(u,v).  Let  any 
positive  number  e  be  given.  Then  f(u,v)  being  regular  in  closed  regions 
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one lotrliv^,  the  nti;lp  lUorc  exists  a  &  >  0  anti  a  y  In  -I  y  <  J 

such  that.  |f(u,v)  -  f(u,y)|  <■  e  is  true  for  every  v  on  a  particular  con¬ 
tour  for  which  0  <"  |v  -  y|  <r  6  ,  Also,  f(u,v)  being  a  regular  function  of 
V  in  regions  enclosing  the  trips  (-1,1),  tVic  oscillations  of  f(u,y)  with 
respect  to  y  are  finite,  and 


f(u,y) 


is  true  for  every  y  in  -1  <  y  <  1.  Thus 
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But  the  argument  on  the  right  hand  side  of  this  equation  is  just 

f  Egy)  dy 

J  ^ 

-1 

V7here  E|(y)  is  given  on  the  right  side  of  (4)  --  this  being  an  error  in 
our  Gaussian  integration  which  we  assume  to  be  able  to  estimate  (a  method 
of  obtaining  a  bound  on  this  error  will  be  illustrated  later)  and  make 
appropriately  small.  Under  these  assumptions  we  can  write  (5)  as 
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since  E^  -  Eg  will  then  only  be  a  negligible  second  order  error. 


Combining  the  results  of  the  above  equations,  we  have 
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where  the  x  are  zeros  of  the  Legendre  polynomial  P  (x^)  and  the 

J  i  J  i 

are  the  corresponding  weights.  Under  the  assumption  that  v/e  can  appropriately 

bound  each  of  the  E^  in  (7b)  we  have  neglected  second  and  higher  order 

errors  on  the  right  side  of  (7a). 


The  practical  use  of  the  error  term  in  (7b)  depends  on  the  fact 
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account  a  finite  number  of  poles  of  f(x  ,  ..,,u  )  lying  in  the. 
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complex  domain  of  u  ,  and  we  show  this  briefly  at  the  end  of  this  section.). 
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It  is  convenient  to  chose  as  contour  a  circle  of  sufficiently  large  radius 
R  and  to  employ  the  asymptotic  value,  |u|  -+00  of  Q  (u)/!-*  (u); 


m 


/r.v  2^(ml  -2m-l  r  2m^+5m^-m-l  -2  /  -4v-, 

=  lisrj'fatryr  “  fi "  ■ 


As  a  simple  illustration  consider 
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We  find  on  taking  only  the  first  term  of  the  asymptotic  expansion  for 
Q  in)/?  (u)  that  the  modulus  in  E,  is  dominated  by 

4y^K^  exp(-[(2m-4)logR“|y|R]3  <  4K^  exp(-[(2m-4)logR-R] } 

on  a  large  circle  of  radius  R,  the  factorial  constant  being  denoted  by 

K  .  The  least  value  of  the  function  in  brackets  qua  function  of  R  is 
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to  obtain 


e.g.  if 


lE^I  <  4K^  exp[-[(2n-2)(log(2n-2)  -  l)]) 


n  =  5,  |E  1  <  10.1  X  lo"’^ 


n  =  6,  |E  1  <  3.18  X  10 
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Hence  if  we  set 
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we  have  for 
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m  =  7,  n  =  6,  iR  I  <  1.13  X  lO"^  . 


These  error  bounds  could  be  improved  by  a  more  careful  analysis  but  they 
are  simply  obtained  and  not  unduly  pessimistic,  the  actual  errors  being 
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If  f(u,y)  is  a  meromorphic  function,  the  contribution  of  the 
residues  at  the  poles  of  f(u,y)  to  the  right  side  of  (2)  may  be  supposed 
to  be  8(x,y)  and  (4)  is  then  replaced  by 
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The  right  side  of  (9)  vanishes  and  the  integration  formula  is 
exact  if  f(u,y)  is  a  meromorphic  function  (of  u)  and  if  for  every  y  in 
-1  <  y  <  1 

|f(u,y)|  <  0(lul^“^),  |ul  -►  00 


uniformly  with  respect  to  the  argument  of  u.  For  example  let  f(u,y)  be 
meromorphic,  satisfying  the  above  inequality.  Let  f(u,y)  have  simple  poles 
(not  located  on  the  strip  (-1  <  u  <  1))  a^(y)  with  residues  A^(y).  We 

then  have 


m 


(10) 


r  f(x,y)  dx  =  ^  aj^£(x^,y)  '  2  V 
'-1  1=1  k 


\(y)  Q^[°T,(y)l 


From  the  expansion  (8)  for  that  if  the  small¬ 

est  of  |a^(y)|  in  (lO)  is  sufficiently  large,  then  by  increasing  m  by  1 

r  1  1  ^ 

we  would  reduce  the  sum  on  the  right  by  a  factor  |2j^(yy| j  ’ 

increase  the  accuracy  in  our  numerical  integration  by  increasing  the  number 

of  evaluation  points. 


If,  however,  some  of  the  poles  close  to  the  strip 

-1  <  Ri  u  <  1,  we  may  need  a  large  number  of  points  m  to  reduce  the 
effect  of  these  poles.  Although  equation  (9)  indicates  an  alternative 
method  of  evaluating  an  Integral  whose  integrand  has  poles  in  the  complex 
plane,  the  residue  g(x,y)  will  in  most  cases  be  a  more  complicated  func¬ 
tion  than  f(x,y)  and  It  will  in  general,  be  better  to  evaluate  the  orig¬ 
inal  Integral  using  a  larger  number  of  points.  This  is  particularly  so  in  n 
dimensions  since  here  the  polos  and  residues  at  the  poles  may  be  functions  in 


n-1  variables. 
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for  A, 
k 

changes 

cussion 


„  .  f  I  i-l 

liy  writing  ,x 

g(x\  .  .  .  ,x^,  .  .  .  for  g( 

corresponding  1;o  those  between 

of  meromorphic  functions  can  o 


^ .  jX*^')  for  a|^(y)  '-WkI  slnii iarly 
x,y),  and  luaking  the  other  appropriate 
equations  (6)  and  (7)>  the  above  dis- 
asily  be  extended  to  n  dlmenslonr. . 


2)  Gauss -Laguerre  and  Gauss-Hermite  formulas 

Denoting  the  zeros  of  the  Laguerre  polynoraials  L  (x'^)  by  xp 

m .  k  . 

.  J  J 

and  the  corresponding  weights  by  c^‘|^  we  illustrate  in  what  follows  a 

method  of  obtaining  a  bound  on  the  error 
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by  employing  contour  integration  and  asymptotic  expansions.  The  method 
used  in  the  previous  section  is  applicable  here,  but  some  modification  is 
needed  since  the  range  of  integration  is  infinite  in  the  Gauss-Laguerre 
and  Gauss-Hermite  formulas.  Vie  consider  the  analysis  of  only  the  two  dimen- 
sinal  Laguerre  formulas  in  detail,  from  which  we  can  then  write  dowi  the 
corresponding  result  in  n  dimensions.  The  Gauss-Hermite  analysis  is 
similar  and  may  be  omitted. 


Assume  that  the  integral 

I  P  e  f(x,y)  dx  dy 

J 

0  0 

exists  and  that  its  value  is  independent  of  the  order  of  integration.  i\ssume 
further  that  exp(-x-y)  f(x,y)  -»-oo  if  either  x,  or  y,  or  both  become 


infinite,  while  the  other  remains  zero  or  positive,  and  that  f(x,y)  does 
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not  have  any  singularities  in  any  finite  region  of  integratioiiio  Hence 
given  any  e  >  0  (Here  €  may  be  taken  to  be  a  negligible  fraction  of 
our  allowable  error  in  the  computation  of  the  value  of  the  integral.) 
there  exists  an  depending  only  on  e  such  that 


(12) 


r  £(x,y)  dx 


<  € 


for  each  y  >  0  and  for  each 
chosen  depending  only  on  e, 


Ml  >  11^2^ •  A  value  of  may  similarly  be 

for  integration  with  respect  to  y. 


Let  the  zeros  of  the  Laguerre  polynomial 
them  lie  within  the  strip  (0,Mj^)  of  the  real  u  axis.  Let  C.^  be  a 
contour  enclosing  the  strip  (0,M^)  such  that  u  on  satisfies 

|u|^^  <  <  |u|  -  6,  where  0  <  6  <  ^  and  k  >  1  are  otherwise, 

arbitrary  positive  numbers.  Then  if  f(u,v)  is  a  regular  function  of  u 
in  a  region  enclosing  this  contour  for  every  complex  number  v  in  a  simi¬ 
lar  region  of  regularity,  we  can  write, 

^  f(x^,y)  L^(x)  ^  p  f(u,y)  L^^(x) 

ial  C, 


where  y  is  in  0  <  y  <  Mg  >  Mg^.  On  multiplying  this  equation  by 

e  ,  integrating  with  respect  to  x  over  and  Interchanging  the 

order  of  integration  in  the  double  integral  we  obtain 
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We  now  proceed  to  find  the  asymptotic  expansion  of  the  inner 
integral  on  the  right  of  this  equation.  Our  asymptotic  sequence  will  be 
(^),  |u|  -►oo;  hence  any  term  which  is  0[|u|^  exp(- |u|  ]  where  r  and 

k  (k  >  l)  are  arbitrary  positive  numbers,  will  not  contribute,  towards  our 
asymptotic  expansion  and  is  "asymptotically  equal  to  zero"  with  respect  to 
our  asymptotic  sequence.  This  can  be  written  more  concisely  as  0(juj^  exp 
(-|u|^A))%o;  (i)  . 


Our  first  step  is  to  show  that 
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If  the  point  u  is  not  on  the  axis  of  Integration, 
For  suppose  u  =  or  +  iw,  where  a  and  to  are  real 
r  being  arbitrary  positive  numbers.  Then 


this  can  easily  be  shown. 
»  ^  a  and 


L  (x)  dx 
m  ' 


Now  for  k  >  0  but  otherwise  an  arbitrary  positive  number,  it  is  readily 
shown  that 


and  therefore 
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If,  on  the  other  hand  u  is  on  the  positive  x-axis,  we  evaluate 


the  Cauchy  principal  value  of  the  integral: 
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where  w  is  some  number  in  0<5<'rM,  a  and  r  being 

I  u  I  “  "  “ 

arbitrary  positive  numbers.  From  the  above  analysis  we  see  that  the  first 
and  last  integrals  on  the  right  are  asymptotically  equal  to  zero  with  reS" 
pect  to  our  asymptotic  sequence,  (— ) . 


To  show  that  the  center  integral  is  asymptotically  equal  to  zero 
with  respect  to  our  asymptotic  sequence,  we  let  x  =  u  t  when  x  >  ti, 

and  X  =  u  -  t  when  x  <  u  in  the  integral  to  obtain 
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The  integral  on  the  right  obviously  exists  since  the  integrand  is  bounded 
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mean  value  theorem  of  the  differential  calculus,  given  any  t  in  0  <  t  <  w 
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Left  the  msHimum  value  of  this  expression  in  «  w  <  A’  <  w  be  at  the  point 
x'  s  u  +  A',  Honco  on  replaeing  the  above  integral  on  the  right  by  the 
maxlmuffl  value  of  the  integrand  in  the  interval  times  w>  we  find  that 
(15)  is  satisfied  even  if  the  point  u  is  on  the  positive  real  axiSe 
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whenever  <M^<|u|-6,  k>l,  |u|— >-oo,  tbe  Inl 

gral  on  the  right  being  aasnmed  to^a  Cauchy  principal  value  when  tbe  point 
u  lies  on  the  positive  real  axis. 

Expanding  the  denominator  in  the  above  integrand  in  powers  of  x/u 
we  note  that 
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the  first  m  terms  of  our  expansion  Integrating  to  zero  due  to  the  ortho¬ 
gonality  property  of  the  Laguarre  polynomials.  Thus  under  the  conditions 

stated  above 
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By  equation  (12)  and  the  above  arguments  we  can  also  neglect  tlio 


contribution  of  the  integration  from  to  co  of  the  left  hand  side  C'f 

(15)  to  our  asymptotic  expansion.  Thus,  with 
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|u  I  -►00  .  Assuming  we  can  suitably  bound  each  on  the  right  side  of 

(19b)  we  have  neglected  second  and  higher  order  errors  on  the  right  side  of 

Cl9a). 


We  now  return  again  to  our  two-dimensional  model  for  simplicity. 

If  f(u,y)  is  a  meromorphlc  function  with  poles  in  the  complex  u- plane,  the 
contribution  of  the  residues  at  the  poles  of  f(u,y)  to  our  contour  inte¬ 
gration  may  be  supposed  to  be  g(x,y)i  and  (16)  Is  then  replaced  by 

(20)  I"  e”^  f(x,y)  dx  -  ^  cj  £(x^,y)  -  T  e"*^  g(x,y)  dx 
0  lal  0 

while  0  <  y  <  . 

m  m  2 

The  right  side  of  (20)  vanishes  and  the  integration  formula  i© 
exact  if  f(u,y)  is  a  maromorphlc  function  of  u  and  if  for  every  y  in 
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uniformly  with  respect  to  the  argument  of  u.  For  example,  let  f(u,y)  be 
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meromorphic,  satisfying  the  above  inequality, 
(not  located  on  the  strip  0  <  u  <  oo)  a,  (y) 

“  K 

have 


Let  f(u,y)  have  simple  poles 
with  residues  A^(y).  We  then 
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/  e  ^  f(x,y) 
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dx 


If  the  a^(y)  sufficiently  far  from  the  origin,  the  contribu¬ 

tion  of  the  residues  to  the  integrand  will  be  negligible.  Since,  as  m  ->oo, 
the  right  hand  side  of  (2l)  becomes  unbounded*,  we  may  want  to  use  formula 
(20 )  to  evaluate  the  original  integral  to  the  accuracy  desired.  This,  however, 
may  become  a  very  complicated  problem,  since  the  function  g(x,y)  will  in 
general  be  more  complicated  than  the  function  f(x,y). 

The  above  discussion  applies  equally  to  the  Hermite  formulas.  It 
can  a  Iso  easily  be  extended  to  the  n-dimensional  case.  In  the  n-dimensioaal 
case  the  poles  and  residues  at  the  poles  may  be  functions  in  (n^l)  variables, 
as  we  have  already  noted  in  our  treatment  of  the  Gauss -Legendre  formulas. 


In  closing,  we  note  that  by  our  procedure  we  can  obtain  error  bounds 
for  multiple  integrals  over  combinations  of  regions,  e.g,  for  integrals  of  the 
form 
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by  a  repeated  Gaussian  Integration  formula*. 


By  our  previously  developed  formulas, 

_  _  2 

p  poo  poo  lOOe  J  (Jy)cos  z 


C  0  “00 


u  -  4u  +  lo4 


■■ 


+  . . . K  du 


|u|. 


>00 


^  ^  I 

2  2x1  J  J  J 

Cg  -«  -1 


00  pi  lOOe  Jq(-|v)cos  2 


-  4x  +  104 


■’*  ‘**1^  *  •••} 


dv 


V  -*■00 


£  <v 

3  2x1 


ii  /  /  / 


00  pi  lOOe  ^J^(^y)cos  w 


C,  0  -1 

5 


O' _ 

x^  -  4x  +  104 


dx  dy 


/  X  nl 

[  n  2n+l 
2  w 


dw 


w  -►00  . 


By  taking  u  sufficiently  large  we  observe  that  =  0.  The  de- 
nominator  of  tha  Integrand  as  a  function  of  u  has  poles  at  the  points 
and  ttg  (ttj^  »  2  +  101,  ag  »  2  -  lOi),  where  ja^l  >  10.  Both 
cos  z  are  less  than  or  equal  to  one  In  magnitude.  Hence  expanding  the  sum 
In  the  residues  of  the  right  hand  side  of  (7)  asymptotically,  l.e.  using  the 
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asymptotic  expansion  for  Q,(a  )/P.(a^),  we  find  that  the  dominating  term  of 
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the  error  is  less  than 
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=  5.14  X  10'^  for  j  =  2. 


O', 


For  X  in  -1  <  x  <  1,  100/(x^  -  4x  +  104)  <  1.  ^qC^v)  has 
maximum  value  for  v  imaginary.  Along  an  imaginary  axis 


l•Jo(4'')|  =  '=°sh  (1  +  o(]^))  • 


Hence  considering  only  the  first  term  of  the  asymptotic  expansion  indicated 
in  Eg  above  together  with  the  dominating  term  of  the  asymptotic  expansion 
for  |jQ(iv)l,  V  imaginary,  we  find  that  the  dominating  term  in  Eg  has  a 
minimum  value  on  a  large  circle  jv|  =  4m  -f  1.  Here 


=  9*^5  X  10  ^  with  m  =  4 


Similar ly 

iE  I  <  ^  (n  +  1)  •  =  1.56  X  10  ^  with  n  =  5» 


The  overall  estimate  of  the  error  in  the  triple  sum  is  jE^j  +  jEgj 

+  |E^1  <  1.4  X  10  .  The  actual  error  when  evaluating  the  integral  by  a 

-6 

numerical  integration  formula  is  7*8  x  10  ,  indicating  that  our  error  bounds 

are  not  too  pessimistic;  we  could  have  improved  them  by  estimating  our  inte¬ 
grals  more  carefully. 


Notice  that  in  our  last  example  we  have  used  our  formulas  for  error 
bounds  to  find  the  required  number  of  points  for  achieving  a  desired  accuracy 
of  numerical  integration. 
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CHAPTER  VI 

SUMMARY,  AND  SOME  UNSOLVED  PROBLEMS 

This  Chapter  states  unsolved  problems  and  summarizes  the  results 
presented  in  the  thesis. 

If  we  survey  in  retrospect  the  development  of  numerical  integration 
formulas  we  discern  a  rapid  advancement  in  the  one -dimensional  cases  during 
the  time  of  Gauss  and  Chebychev.  With  the  help  of  orthogonal  polynomials  the 
results  of  these  authors  have  been  extended  to  many  different  ranges  of  inte¬ 
gration  for  many  different  weight  functions.  Indeed,  every  set  of  orthogonal 
polynomials  provides  a  new  integration  formula. 

For  a  number  of  reasons,  numerical  integration  formulas  in  higher 
dimensions  have  developed  more  slowly.  Before  the  development  of  electronic 
digital  computers  (around  19^7)  no  practical  use  could  be  found  for  an  inte¬ 
gration  formula  in  higher  dimensions.  The  availability  of  powerful  computing 
devices  has  however  made  possible  the  study  of  functions  of  several  independ¬ 
ent  variables,  and  new  attempts  at  obtaining  numerical  integration  formulas 
in  higher  dimensions  are  being  made. 

Among  the  formulas  available  in  higher  dimensions  are  those  over 
rectangular  regions  that  can  be  obtained  from  one-dimensional  integration 
formulas  by  taking  products  of  integrals  over  each  separate  variable.  As 
the  number  of  dimensions  increases  the  number  of  evaluation  points  required 
by  these  formulas  Increases  so  rapidly  that  the  economic  limit  on  even  the 
most  powerful  digital  computers  is  very  quickly  reached. 

An  estimation  of  the  minimum  number  of  points  required  by  an  inte¬ 
gration  formula  of  precision  p  tells  us  that  the  integration  formulas 
resulting  from  taking  products  of  Integrals  over  each  separate  variable,  are 
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extravagant  in  points  when  the  number  of  dimensions  n  is  large;  that  is, 
integration  formulas  exist  which  require  less  than  points  to  inte¬ 

grate  an  arbitrary  monomial  of  degree  <  p  over  an  arbitrary  region  in  n- 

B 

space. 


This  class  of  "minimum-point”  formulas  is,  however  not  easily 
obtained.  We  have  seen  earlier  that  the  main  difficulty  is  the  solution  of 
systems  of  non-linear  algebraic  equations.  In  the  one -dimensional  case 
there  exists  a  beautiful  theory  of  orthogonal  polynomials  which  enables  one 
to  obtain  numerical  integration  formulas  with  little  difficulty.  This  theory 
is  readily  extended  to  rectangular  regions  in  n-space,  where  it  enables  one 
to  obtain  -  point  formulas.  Unfortunately  no  such  theory  exists  as 

yet  for  the  class  of  minimum-point  formulas.  The  development  of  such  a  theory, 
or  else  the  development  of  effective  methods  of  solving  a  system  of  non¬ 
linear  algebraic  equations  would  solve  the  problem  of  obtaining  minimum-point 
numerical  integration  formulas  of  precision  p  over  regions  in  n-space. 

A  great  deal  of  work  on  obtaining  numerical  integration  formulas  in 
higher  dimensions  has  been  done  by  P.  C.  Hammer  and  A.  H.  Stroud.  They  have 
given  numerical  Integration  formulas  up  to  precision  3  £0^  various  regions 
in  n-space,  and  formulas  up  to  precision  5  £01^  n-cube  and  the  n-sphere. 

In  cases  where  a  transformation  is  available  between  two  regions  they  have 
given  a  method  for  extending  numerical  integration  formulas  from  one  region 
to  the  other  region.  They  have  introduced  and  applied  the  concept  of  fully 
symmetric  regions  which  reduces  considerably  the  number  of  simultaneous 
non- linear  algebraic  equations  we  need  to  solve  to  obtain  numerical  Inte¬ 
gration  formulas  over  regions  In  n-space. 


Numerical  Integration  formulae  of  precision  3,  5  «nd  9  for  arbit 
rary  fully  symmetric  regions  in  n-space  are  developed  In  this  thesis.  The 
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Apart  from  conceivable  fvally  symmetric  regions  and  weight  functions  for  which 
equations  such  as  (lO)  on  page  UU  do  not  hold  the  method  by  which  these  formulas 
were  obtained  will  lead  to  fully  symmetric  numerical  integration  formulas  of 
precision  4k  +  1.  The  non-linear  unknowns  in  the  simultaneous  non-linear  alge¬ 
braic  equations  set  up  to  obtain  numerical  integration  formulas  are  shown  to  be 
the  zeros  of  polynomials  in  one  variable  orthogonal  over  the  fully  symmetric 
region  with  respect  to  the.  fully  symmetric  weight  function. 


To  overcome  the  restriction  the  above  equation  imposes  on  p 
(p  ^  4k  -  1)  we  turn  to  the  repeated  Gaussian  formulas j  noting  that  for  moder¬ 
ate  n  and  large  p  these  formulas  are  not  unduly  extravagant  In  evaluation 
points  inasmuch  as  they  produce  greater  accuracy  for  large  p.  A  transforma¬ 
tion  is  given,  by  which  we  can  transform  the  integral 
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into  an  integral  over  the  n-cube.  This  transformation  serves  to  broaden  the 
range  of  application  of  the  Gaussian  formulas,  and  suggests  that  instead  of 
looking  for  numerical  integration  formulas  over  arbitrary  (bounded)  regions 
in  n-space  we  may  restrict  ourselves  to  the  simpler  problem  of  looking  for 
numerical  integration  formulas  for  the  n-cube,  varying  only  the  weight  func¬ 
tion.  For  completeness,  a  2(^“)'^-point  numerical  integration  formula  of 
arbitrary  high  precision  p  is  developed  for  the  finite  and  infinite  n- 
sphere,  where  n  >  2. 

Error  bounds  are  necessary.  The  majority  of  the  a  priori  estimates 
require  an  estimate  in  the  form  of  a  (p+l)  -th  derivative  of  the  integrand 
while  in  most  cases  empirical  estimates  are  more  easily  obtained  than  moder¬ 
ate  or  high  order  derivative  estimates.  In  this  thesis  we  have  utilized  con- 
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tour  integration  and  asymptotics  to  obtain  error  bounds  for  repeated  Gaussian 
integration  in  the  case  where  the  integrand  is  a  meromorphic  function  of  n 
complex  variables. 

The  Appendix  of  the  thesis  consits  of  (A)  a  number  of  examples  of 
numerical  integrations  over  regions  in  l|.-apace  and  (b)  Table  V  which  is  a 
tabulation  of  the  zeros  and  weights  of  four  particular  fully  symmetric  numer¬ 
ical  integration  formulas  of  precision  9. 

Computers  may  ask  what  positive  recommendation  can  be  made.  The 
field  has  not  yet  been  explored  to  the  extent  that  precise  prescriptions  are 
possible,  but  it  may  be  useful  to  give  some  broad  indications. 

In  the  one-dimensional  case  it  is  easy  to  find  examples  of  inte¬ 
grands  favorable  to  one  type  of  numerical  integration  formula  and  unfavorable 
to  another.  This  suggests  that  since  we  are  unlikely  to  find  omnibus  inte¬ 
gration  formulas  in  the  one  dimensional  case,  we  are  even  less  likely  to  find 
such  formulas  in  n  dimensions.  In  practice  we  use  Gaussian  integration  in 
the  one  dimensional  case,  tolerating  a  certain  amount  of  extravagance  if 
necessary.  Extravagance  is  unfortunately  less  tolerable  in  n-space,  and  in 
order  to  keep  the  number  of  points  at  a  feasible  level  the  choice  of  the 
integration  formula  will  depend  to  a  large  extent  on  the  integration  problem. 

Higher  dimensional  integrals  often  come  up  in  computing  centers-- 
mainly  from  nuclear  physics  and  chemical  engineering.  No  attempt  has  been 
made  to  collect  examples  of  these  or  to  classify  them,  although  this  would 
be  desirable. 

At  present,  moderate  precision  in  higher  dimensional  Integration 
is  relatively  easy  to  obtain.  High  precision  is  also  relatively  easy  to 
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obtain,  but  may  sometimes  be  costly.  This  depends  largely  on  the  C07ii[:)lexity 
of  the  integrand,  the  number  of  dimensions,  and  the  precision  desired.  The 
error  of  numerical  integration  is  sometimes  tractable,  though  usually  it  is 
intractable  by  reason  of  the  complexity  of  the  integrand.  Usually  the 
numerical  integration  is  repeated  with  a  formula  of  higher  precision  and 
the  results  of  the  two  integrations  are  compared.  This  procedure  is  some¬ 
what  unsatisfactory  since  it  can  be  misleading. 
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APPENDIX  A 

EXAMPLES  OF  NUMERICAL  INTEGRATIONS 


In  this  Appendix  we  give  examples  to  illustrate  the  accuracy  of 
some  of  the  numerical  Integration  formulas  developed  in  this  thesis.  The 
integrals  were  evaluated  on  the  I.B.M.  1620  at  the  University  of  Alberta. 

-1  -1  -1  -1 

Tb§  results  of  this  example  are  given  in  Table  VII,  This  Table 

contains} 

(1)  the  value  of  a 

(ll)  the  exact  value  of  the  integral 

(ill)  relative  errors  E(p3),  E(p5)»  E(p9A),  E(p9B),  By  relative 

error  is  meant  the  absolute  value  of  the  actual  error, 

8-1  divided  by  I  ,  where  S  is  obtained  by  numerical 
m  a  a  m 

Integration. 

The  notation  1(p3) . . .E(p9S)  is  explained  below, 

E(p3)  is  the  relative  error  of  the  8-polnt  formula  of  precision 
3  for  the  4-§ube  given  on  pages  4l-h3. 

E(p3)  Is  the  relative  error  of  the  33-point  formula  of  precision 
3  for  the  4-eub@  given  on  pages  4l-43> 

E(p9A)  end  E(p9B)  are  relative  errors  of  the  177-polnt  formulas 
of  precision  9  given  on  pages  42-46  and  tabulated  on  page  109  Tablir^  V 
(Appendix  B).  Formula  A  Is  the  formula  tabulated  on  the  left  hand  side  of 
this  tablet  formula  B  Is  tabulated  on  the  right  hand  side. 


I  SI 


[w  +  X  +  y  +  2  +5]°^  dw  dx  dy  dz 


aKorTM03T«i  oADiisamm  ito  saawAxa 

*  ^ 


i 


39  X^tittsotiZ  9rfJ  o3  esiq/ufix#  svig  w  xibn&qqA  sl/fj  «I 

.aiadrfa  eld3  al:  bsqolsvab  aeiiifnaoJ  noJ:Sa’ig®:ifll  J8ol"isau/«  »dJ  3o  omoe 
io  \faia:tavlnU  arts  Ja  OSdl  .M.a.I  arfi  no  baaaulsva  »i®w  ala-xgs^ni 


u 


r?<nT  elrfT  ,IIV  ^Irfa'X  ni  navig  alqmax^  airfg  |o  a^Tityid?  e^fX 


«  to  »wjt«v  HHS  (i) 
9tt3  3o  ooi^iv  ^aaJi»  #ff3  ,  X  (xx) 

ivf3u«  .(«t!9)a  .(At’<i)a  .(2,)B  ,(jq)s  *wt»  {in) 

j'jO'tTca  lisuiaa  adi  to  awlAv  »ii>Xoida  »41  Jt/iMin  ti  lO'S'sa 
'i'i&mrf  hani»3<So  ai  8  j  bobXvXfe  I  fi 

W  JD  JO  fli 


•  walad  bwXA.rqjf«  at  (a<^<l)a.  ^(5:q)8  flol3fiion  erfT 
olixa^^jq  to  aXi>fr»tot  anlsq-S  sd3  to  rorta  oyXiali?  orfs  iX  (Cq)I 

e»s«q  00  tmv^$  »dvo*4  trfa  ^fot  5 

to  ali/mol  So  avX'Jsloit  »fis  bX  (l!q)®  ® 

Mi  . 

P^lig  «e  otvXa  »^t/o»4  »d3  *ot  5 

■jij 

t-aol  5nXoq-'^p  ida  So  «Troito  $vi3»I»%  «-i8  boa  (^q)2 

V  A-idat  So  (?0I  «8«q  oo  bffSiiXtfdBjJ  km  $4*®^  m$»q  00  w#¥lt  ?  floX«Xo»i[q  to 

t5  , ®^'' 

'  ^  -  b^'OaXwda^  aXw«ot  «.djr' iltmol  . 

/'’■ivv;- '  -abX®  boad  3dfflX:*  ads  no  bajtfili/daJ  si'  Ali/myni  lAlda® 


100 


TABLE  VII  [EXAMPLE  l)] 


a 

I 

(X 

E(p3) 

e(p5) 

e(p9A) 

e(p9b) 

-.5 

7.3166065 

1.6x10  ^ 

,  -4 
2.2x10 

7.9x10"^ 

l.Oxlo'"'' 

.5 

35.525901 

2.4x10'^ 

3.7x10'^ 

1.2xl0'^ 

1 

c 

1.5 

182.49692 

-4 

1.2x10 

3.3x10' 

4.ixio'^ 

5.1x10'^ 

2.0 

421.33189 

3.6x10"^  * 

3.0xio'^  * 

~6 

7.0x10  * 

“6 

2.3x10  * 

2.5 

983.59047 

l.8xio"^ 

6.9x10“^ 

5.0xio'^ 

2.6x10'^ 

4.0 

13,276.768 

3.6x10'^ 

1.9x10"^  * 

2.3x10"^  * 

-6 

2.1x10  * 

6.0 

479,213.79 

3.9xlo’^ 

4.4x10"^ 

2.9x10“^  * 

4.0xio'^  * 

8.0 

19,236,265 

1.2x10"^ 

8.0x10'^ 

2.1xl0’^  * 

-6 

1.0x10  * 

10.0 

840,845,570 

2.5xl0"^ 

3.7x10'^ 

3.2xl0"^ 

l.Oxlo'^ 

2)  I  =  /  /  r  /  _ _ - 

Ct-  J  J  J  J  S^)®'  ' 

4-sphere 

The  results  of  this  example  are  given  in  Table  VIII,  a,  I  » 

a 

E(p5)>  E(p5),  E(p9A)  and  E(p9B)  are  defined  as  in  the  previous  example. 
The  precision  3  5  formulas  are  located  on  pages  41-43.  The  precision 

9  formulas  are  given  on  pages  42-46  and  tabulated  on  page  II9  of  Table  V 
(Appendix  B), 

ERl"^,  ER2'^  and  ER3’^  are  relative  errors  of  the  2(m)’^-polnt 
formulas  of  precision  2m- 1  over  the  n-sphere  derived  in  section  3)  of 
Chapter  IV. 


* 


Theoretically,  the  relative  errors  should  have  vanished  in  these 
eases.  The  error®  given  are  due  to  an  accumulation  of  round-off 
errors;  8  significant  figures  were  carried  In  the  calculations. 
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See  the  footnote  on  page  100. 
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ERl  Is  the  relative  error  of  the  2-point  formula  of  precision  1; 

ER2'*’  is  the  relative  error  of  the  32“Point  formula  of  precision  3* 

ER3^  is  tVie  relative  error  of  the  l62-polnt  formula  of  precision  5- 

The  two  examples  given  so  far  serve  to  indicate  that  branch  points 
in  the  Integrand  have  a  less  destructive  effect  on  the  accuracy  of  numerical 
integration  than  poles  in  the  complex  plane. 

Note  also  that  formula  A  of  precision  9  tends  to  give  slightly 
better  results  than  formula  B.  The  most  likely  reason  for  this  is  that  the 
weights  in  formula  A  have  a  better  distribution  of  magnitudes  than  those 
in  formula  B.  This  suggests  that  we  could  achieve  better  accuracy  in  our 
numerical  integration  by  first  carrying  out  the  sums  that  are  multiplied  by 
the  smallest  weights.  The  accuracy  of  the  2(m)^  -point  formulas  are  sur¬ 
prising.  The  extremely  large  number  of  points  required  would  however,  make 
them  useless  when  the  number  of  dimensions  is  large;  this  rapid  growth  of 
required  evaluation  points  does  not  occur  for  the  other  formulas. 


3)  I  .  r'  r'  r’  r" 

j  j  j  j 


[w4;X+y+z4'5]'^dw  dx  dy  dz 

[ ( 1-w^) ( 1-x^) ( 1-y^) ( 1“Z^) ]2 


The  results  of  this  example  are  given  in  Table  IX.  In  this  Table, 

a,  I  ,  E  (p3),  E(p5),  E(p9A)  and  E(p9b)  are  defined  as  in  Example  l).  The 
a 

precision  3  and  5  formulas  are  located  on  pages  The  precision  9 

formulas  are  given  on  pages  42-46  and  tabulated  on  page  129  of  Table  V 
(Appendix  B).  Formula  A  is  the  formula  tabulated  on  the  left  hand  side  of 
this  Table;  formula  B  is  tabulated  on  the  right  hand  side. 
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TABLE  IX  [EXAMPLE  5)] 


a 

I 

a 

E(p3) 

e(p5) 

e(p9A) 

e(p9b) 

-.5 

45.159877 

4.2xl0"^ 

-4 

1.2x10 

2.1x10*^ 

5.7x10"^ 

.5 

215.45965 

5.2xl0"^ 

“4 

1.0x10 

6.7x10"^ 

1.4x10*^ 

1.5 

1,122.2089 

-4 

2.7x10 

2.2x10*^ 

4.5x10*^ 

4.1x10*^ 

2.0 

2,630.0118 

1.3x10"^  * 

1.1x10*^  * 

1.3x10*^  * 

1.1x10*^  * 

2.5 

6,258.3709 

3.8x10”^ 

1.8x10*^ 

5.6x10*^ 

4.6x10  ^ 

4.0 

91,125.664 

6.9x10-3 

3.7x10"^  * 

4.6x10  ^  * 

4.6x10*^  * 

6.0 

3,739,737.7 

1.7X10”^ 

1.2x10*^ 

1.1x10*^  * 

1.2x10*^  * 

8.0 

173,555,240 

1.9x10"^ 

1.8x10  ^ 

1.4x10*^  * 

1.2x10*^  * 

10.0 

8,836,719,600 

3.5x10*’' 

7.3x10’^ 

1.0x10*^ 

2.4x10*^ 

h) 


1  ^  r°°  r*”  r*^  expC-w^-x^^v^-z^)  dw  dx  dv  dz 

^  J  J  J  J  [  l+w^+x^+y^+z^]^ 


-00  “00 


The  results  of  this  example  are  given  in  Table  X.  All  notations 
are  the  same  as  those  for  Example  l).  The  integration  formulas  of  precision 
9  are  tabulated  on  page  139  of  Table  V  (Appendix  b). 

The  results  of  Example  3)  (Table  IX)  are  better  than  the  results 
of  Example  (Table  X).  This  is  largely  due  to  the  integrand  being  similar 
in  complexity  in  the  two  examples,  while  the  region  of  integration  is  larger 
in  Example  it-).  We  note  again,  that  the  destructive  effect  on  the  accuracy 
of  the  numerical  integration  is  much  more  due  to  poles  in  the  complex  regions 
of  the  Integrand  than  to  branch  points.  Also,  formula  A  of  precision  9 


* 


See  footnote  page  100. 
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gives  somewhat  better  results  than  formula  B;  this  being  largely  due  to 
the  weights  In  formula  A  having  a  much  better  distribution  of  magnitudes. 


TABLE  X  [EXAMPLE  4)] 


a 

1 

. -  . . a . 

B(PJ) 

e(p5) 

e(p9a) 

e(p9B) 

“3*4 

1,092.7610 

5.7X10“^ 

2.3x10*^ 

2.4x10“^ 

o 

1 

403,61046 

4.5X10“^ 

1.2x10“'^' 

-A 

4.1x10  * 

4.8x10“’^  * 

-2.5 

223. 01810 

3.1X10“^ 

4.5x10“^ 

1.8x10“^ 

1.1x10“^ 

-1*5 

55.295497 

7.3X10"^ 

l.lxlO 

«4 

1.5x10  ^ 

1.2xlo“^ 

“1.0 

29. 600802 

l.OxlO^'^  * 

3.4x10‘’'^  * 

2,0x10*“*^  * 

2.4x10“*^  * 

•  .5 

16.67^363 

2.55^10“^ 

1.3x10 

2.4xl0“^ 

2  7xl0"^ 

1.5 

2.7013^33 

3.0x10“^ 

“1 

5.2x10“^ 

l.lxlo“^ 

-2 

1.9x10'’^ 

6.3x10“^ 

2.5 

1.3Sa6TT2 

5 . 4x10 

1.5 

2.0x10 

4" 


This  is  due  to  Js  Oadwell  [2^3 1  iu  fehli  fuper  an  Algol 

prOft'Kifl  Is  given  to  evaluate  an  Intugfal  of  the  type  (S)*  Ghapter  IV  (page 
64)  using  llfflpson*®  ful§**  with  rtpeated  blseotlon  to  attain  the  required 
aecuraoy.  In  this  paper  the  author  has,...  In  ©ffeot,  utilised  the  transfornta' 
tlon  given  in  Ghapter  W  of  this  thesis «  We  present  the  results  of  this 
eisample  in  fable  11. 


*  let  the  footnote  on  page  100. 

« jc^+ikh  I  _ 

**  ilmpson^s  lules  /  "  =  |[f^+4£|+ifg+4f^+. . 

0 

where  f^  s  f(jt^+Jh)  . 
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^(p3)>  1'(p5)j  E(p9A),  E(p9B),  ERI,  ER2  and  ER3  have  the  same 
meaning  as  in  Example  2).  ER(S1),  ER(S2)  and  ER(S3)  are  the  relative 

errors  by  the  method  in  [23]  •  i'he  number  of  evaluation  points  is  also  given. 


TABLE  XI  [EXAMPLE  5)] 


Exact  Value 

2.6079550 

Relative 

Errors 

No.  of  Eval¬ 
uation  points 

E(P3) 

2.8x10“^ 

8 

e(p5) 

U.Sxio  ^ 

33 

e(p9A) 

1)  .  6x10  ^ 

177 

e(p9b) 

2.5x10'^ 

177 

ERl 

2.2x10“^ 

2 

ER2 

3.5x10“^ 

32 

ER3 

l.Oxlo"^ 

162 

ER(Sl) 

6. 2x10 

625 

ER(S2) 

6.2xl0“^ 

3,345 

ER(S3) 

^  ~h 

6.2x10 

69,113 

Wote  again,  the  surprising  accuracy  of  the  2(m)  ”  point  formulas 

constructed  over  the  n-sphere,  and  that  the  l77“Poiiih  formula  A  of  prec¬ 
ision  9  gives  slightly  better  results  than  fornmla  B. 

If  proper  Gaussian  formulas  had  been  used  instead  of  Simpson's 
rule  in  [23]  the  given  accuracy  could  have  been  obtained  with  less  than  ^ 
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of  the  number  of  points  used*.  Moreover,  taking  advantage  of  the  full 
symmetry  (as  was  done  in  [2^]l)  of  the  integrand  we  could  have  reduced  the 
number  of  points  by  an  additional  factor  of  2^  ■-  16  (an  even  number  of 
evaluation  points  is  assumed  in  each  one-dimensional  integration  formula 
from  which  we  obtain  an  integration  formula  over  a  U-dlmensional  cube). 
Nonetheless,  we  note  the  superior  accuracy  of  the  integration  formulas 
constructed  for  the  particular  regions  over  which  they  are  used  as  compared 
with  integration  formulas  constructed  over  rectangular  regions  and  used 
over  curvilinear  regions. 


The  following  procedure  is  Implied  here.  ^ 

First  apply  the  transformation  theorem  given  in  the  section 

of  Chapter  IV.  The  resulting  Jacobian  of  the  transformation 
breaks  up  into  suitable  weight  functions  with  respect  to  which  the 
Gegenbauer  (Ultraspherical)  polynomials  [I3]  are  orthogonal  over 
the  interval  (-1,1).  The  numerical  integration  formulas  (re¬ 
peated  sums)  resulting  from  using  the  zeros  of  these  polynomials 
as  evaluation  points  can  then  be  used  to  evaluate  the  integral. 

We  have  applied  the  method  of  Chapter  V  to  obtain  the  above  esti¬ 
mate  on  the  number  of  points  required. 
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APPENDIX  D  -  TABLE  V 

ZEROS  AND  WEIGHTS  OF  FORMULAS  OF  PRECISION  NINE 


In  this  Table  we  have  tabulated  zeros  U,  V  (U  =  u,  V  =  v)  and 
weights  A,  ii,  C,  D,  E,  F,  H,  I,  J  (G  =  O)  of  the  precision  nine  integration 
formula(s)  given  on  pages  42-46.  The  zeros  and  weights  are  for  the  four 
following  regions  and  weight  functions; 


Region 

N  -  Cube 
N  -  Sphere 

N  -  Cube 

Infinite  N-Space 


Weight  Function 


1 

1 

N 

i/n  [I  - 

issl 


N 


Tabulation  of  the  zeros  and  weights  is  carried  out  for  Increasing 
dimension  N(N  s  n)  starting  from  4. 

Only  nine  out  of  the  twelve  equations  (page  45)  were  required  to 
compute  the  weights.  The  other  three  equations  were  used  as  checks  to  compute 
the  errors  in  the  computations.  These  errors,  ERl,  ER2,  and  ERJ  are  listed 
following  the  listing  of  the  zeros  and  weights.  Equation  (lx)  page  45  was 
used  to  compute  EKl}  equation  (v)  page  45  was  used  to  compute  ER2;  and 
equation  (iv)  page  45  was  used  to  compute  ERJ.  On  comparing  the  magnitude 
of  the  errors  and  weights  we  find  that  the  largest  of  B  or  C  may  be  in 
error  by  2  in  the  19th  significant  figure.  Since,  however,  the  data  were 
obtained  by  truncating  20  significant  figure  results  to  16  significant 
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figures  we  feel  that  the  largest  of  B  or  C  (and  also  A)  is  correct  to 
18  significant  figures  unless--unless  the  offline  printer  has  made  errors. 
When  in  good  running  order  it  is  estimated  that  the  offline  printer  may  make 
one  error  in  printing  100,000  digits.  Since  the  printer  was  in  good  running 
order  when  the  results  were  printed  it  may  have  made  two  errors.  Any  other 
weight  R  will  have  at  least 

18  -  log  saslSjCi 

logfo  ^ 

correct  significant  figures  since  the  results  were  obtained  in  the  order 
given  on  pages  itl|-  and  1+6. 

Computations  were  carried  out  using  20  figure  floating  point  arith¬ 
metic  on  the  I.B.M.  1620  at  the  University  of  Alberta.  The  printer  referred 
to  above  is  the  I.B.M.  I+07. 
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u 
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= 
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= 
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D 
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N  - 

•  CUBE 

W<  X 

)  t  1 

N  = 

7 

u 

— 

.538469 

310105 

683091 

E-00 

U 

— 

.906179 

845938 

663992 

E-00 

V 

= 

•906179 

845938 

663992 

E-00 

V 

= 

.538469 

310105 

683091 

E-00 

A 

= 

-.132655 

593205 

210437 

E  +  04 

A 

= 

.291049 

998350 

767516 

E+04 

B 

= 

.303966 

168610 

599818 

E  +  03 

B 

= 

-.249520 

144267 

863757 

E+02 

C 

= 

-.440669 

328410 

587296 

E  +  02 

C 

= 

-.925736 

154249 

323833 

E+03 

D 

= 

-.432753 

081799 

139034 

E  +  02 

D 

= 

-.132693 

269800 

104851 

E+01 

E 

= 

.345179 

690556 

703995 

E  +  01 

E 

= 

.264150 

272535 

067009 

E+03 

F 

= 

.362880 

000000 

000000 

E  +  01 

F 

= 

.362880 

000000 

000000 

E+Ol 

H 

= 

-.955745 

920713 

617694 

E-00 

H 

= 

-.614851 

161429 

961826 

E+02 

I 

= 

.126515 

102787 

704048 

E+01 

I 

= 

.780807 

965399 

691215 

E-01 

J 

= 

.197549 

036629 

171333 

E-OO 

J 

= 

.895079 

054575 

156330 

E+01 

U 

.538469 

310105 

683091 

N  = 

E-00 

8 

U 

.906179 

845938 

663992 

E-00 

V 

= 

.906179 

845938 

663992 

E-00 

V 

.538469 

310105 

683091 

E-00 

A 

= 

-.444911 

507755 

310953 

E  +  04 

A 

= 

.152724 

542748 

705916 

E+05 

B 

= 

.887972 

868617 

051137 

E+03 

B 

-.516157 

415937 

315217 

E+02 

C 

= 

-.129465 

591135 

495290 

E+03 

C 

= 

-.412026 

931903 

063791 

E+04 

D 

= 

-.106793 

032805 

860454 

E  +  03 

D 

= 

-.390315 

814064 

160296 

E+01 

E 

= 

.127789 

524754 

506331 

E  +  02 

E 

= 

.966401 

721690 

072914 

E+03 

F 

= 

.725760 

000000 

000000 

E+01 

F 

= 
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000000 

000000 

E+01 

H 

= 

-.278035 

176934 

870602 

E+01 

H 

= 

-.178865 

792415 

988894 

E+03 

I 

= 

.202424 

164460 

326476 

E  +  01 

I 

= 
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274463 
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E-00 

J 

= 
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E-00 
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208862 
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E+02 
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U 
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— 
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V 
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E-00 

V 

= 
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310105 
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A 

= 
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B 
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N  -  SPHERE  W( X )  =  1 


N  =  7 


u 

= 

•385270 

382880 

547445 

E-00 

U 

= 

.719884 

295384 

850745 

E-00 

V 

= 

.719884 

295384 

850745 

E-00 

V 

= 

.385270 

382880 

547445 

E-00 

A 

= 

-.853984 

718832 

534907 

E  +  01 
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= 

-.415796 

381192 

531169 

E+01 
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= 

.199595 

992945 

675215 

E  +  Ol 
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= 

-.228900 

991731 

237094 

E-00 
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= 

-.237383 

980846 

572939 

E-00 
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= 

.735510 

261369 

977054 

E-00 
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= 

-.304649 

099939 

537498 

E-00 
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= 

-.240324 

148253 

265722 

E-02 

E 

= 

-.282494 

203698 

695919 

E-03 

E 

= 

.104633 

170821 

562761 

E-01 
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= 

.348988 

395535 

842131 

E-01 

F 

= 

.348988 

395535 

842131 

E-01 
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= 

-.424149 

455766 

792260 

E-03 
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= 

-.630224 

834043 
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E-01 
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= 

.370551 

990877 

550432 
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744336 
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N  -  CUBE  W(X)  =  1/||  [l-(x^)2]i 

i=l 


N  =  10 


u 

= 

.587785 

252292 

473129 

E-00 

U 

= 

.951056 

516295 

153572 

E-00 

V 

= 

.951056 

516295 

153572 

E-00 

V 

.587785 

252292 

473129 

E-00 

A 

= 

-.460684 

409827 

191913 

E+07 

A 

= 

.617321 

530586 

825619 

E+08 

B 

= 

.825343 

702815 

861806 

E  +  06 

B 

= 

.212346 

463679 

287028 

E+05 

C 

= 

-.267329 

307368 

680100 

E  +  06 

C 

= 

-.127310 

196823 

116432 

E+08 

D 

= 

-.747096 

849081 

404085 

E  +  05 

D 

= 

-.770059 

687081 

264990 

E+04 

E 

s 

.403933 

954186 

221506 

E  +  05 

E 

= 

.218468 

421261 

311042 

E+07 

F 

= 

.374592 

189904 

332084 

E  +  04 

F 

= 

.374592 

189904 

332084 

E+04 

H 

= 

-.601174 

903617 

935006 

E  +  04 

H 

-.282424 

237190 

156354 

E+06 

I 

= 

.700496 

489349 

854726 

E+03 

I 

= 

.102201 

060445 

142595 

E+03 

J 

= 

.531611 

705886 

524743 

E  +  03 

J 

= 

.208736 

562886 

467371 

E+05 

N  =  11 


U 

= 

.587785 

252292 

473129 

E-00 

U 

= 

.951056 

516295 

153572 

E-00 

V 

= 

.951056 

516295 

153572 

E-00 

V 

= 

.587785 

252292 

473129 

E-00 

A 

= 

-.195245 

252659 

718396 

E  +  08 

A 

= 

.319141 

897170 

582023 

E+09 

B 

= 

.328524 

774451 

221983 

E+07 

B 

.127518 

863984 

911087 

E+06 

C 

= 

-.121170 

331362 

268179 

E  +  07 

C 

= 

-.596296 

876033 

778936 

E+08 

D 

= 

-.265516 

842006 

719376 

E+06 

D 

=: 

-.286871 

759671 

712203 

E+05 

E 

.172196 

150673 

967711 

E+06 

E 

.917172 

346017 

679764 

E+07 

F 

= 

.117681 

607189 

556238 

E+05 

F 

= 

.117681 

607189 

556238 

E+05 

H 

= 

-.223203 

696267 

932146 

E+05 

H 

-.104858 

225579 

816855 

E+07 

I 

= 

.192559 

029670 

607639 

E+04 

I 

= 

.280939 

838098 

103085 

E+03 

J 

= 

.167596 

293977 

267900 

E+04 

J 

.674619 

812840 

916112 

E+05 

N  =  12 


U 

.587785 

252292 

473129 

E-00 

U 

=: 

.951056 

516295 

153572 

E-00 

V 

.951056 

516295 

153572 

E-00 

V 

= 

.587785 

252292 

473129 

E-00 

A 

=: 

-.797116 

722394 

239093 

E+08 

A 

= 

.155713 

747352 

090707 

E+10 

B 

.127541 

131177 

860105 

E+08 

B 

.629304 

630672 

447670 

E+06 

C 

=; 

-.530405 

384389 

959464 

E+07 

C 

= 

-.265987 

436316 

841196 

E+09 

D 

-.930936 

485531 

961004 

E+06 

D 
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E+06 

E 
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E  +  06 

E 
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E+08 
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=; 
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349674 

E  +  05 

F 
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e+05 
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= 
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337441 

642134 

E+05 

H 

a: 

-.380102 

112865 

400736 

E+07 

I 

.537726 

251551 
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E+04 
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= 

.784532 

027973 
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E+03 
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056931 
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E  +  04 
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= 

.216545 
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294275 

E+06 
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IN  i  1 

N  -  CUBE  W(X)  =  1/  [|  [l‘(x^)2]2 

1=1 


M  =  1^ 


u 

= 

.587785 

252292 

473129 

E-00 

U 

.951056 

516295 

153572 

V 

= 

.951056 

516295 

153572 

E“00 

V 

= 

.587785 

252292 

473129 

A 

= 

-.315488 

893269 

137486 

E  +  09 

A 

= 

.726085 

263039 

446232 

E+IO 

B 

= 

.485232 

326338 

009118 

E  +  08 

B 

= 

.281378 

151064 

193867 

E+07 

C 

= 

-.225368 

842263 

272287 

E  +  08 

C 

= 

-.114241 

843674 

372207 

E+10 

D 

= 

-.322870 

025537 

895362 

E  +  07 

D 

= 

-.371859 

560181 

517797 

E+06 

E 

= 

.279697 

365693 

962813 

E  +  07 

E 

= 

.145639 

008416 

811419 

E+09 

F 

= 

.116147 

090824 

531336 

E  +  06 

F 

= 

.116147 

090824 

531336 

E+06 

H 

= 

-.288075 

747011 

013266 

E  +  06 

H 

= 

-.135334 

280611 

082157 

E+08 

I 

= 

.152038 

515736 

522043 

E+05 

I 

= 

.221821 

205002 

749606 

E+04 

J 

= 

.166219 

994005 

781593 

E  +  05 

J 

= 

.691875 

254629 

062990 

E+06 

N  =  14 


U 

= 

.587785 

252292 

473129 

E-00 

U 

= 

.951056 

516295 

153572 

E-00 

V 

= 

.951056 

516295 

153572 

E-00 

V 

= 

.587785 

252292 

473129 

E-00 

A 

= 

-.121614 

538388 

747660 

E  +  10 

A 

= 

.326409 

304444 

894127 

E+11 

B 

= 

.181549 

841164 

132206 

E+09 

B 

= 

.118401 

878890 

600800 

E+08 

C 

= 

-.933546 

886616 

863696 

E  +  08 

C 

= 

-.476037 

015372 

616986 

E+10 

D 

= 

-.110985 

471711 

431176 

E  +  08 

D 

= 

-.130760 

563604 

280264 

E+07 

E 

.108302 

647351 

971723 

E+08 

E 

= 

.559122 

990700 

814813 

E+09 

F 

= 

.364886 

847270 

174126 

E+06 

F 

= 

.364886 

847270 

174126 

E+06 

H 

= 

-.101148 

919760 

333124 

E+07 

H 

= 

-.475184 

614893 

298275 

E+08 

I 

= 

.434220 

985545 

958021 

E  +  05 

I 

= 

.633519 

880042 

794231 

E+04 

J 

= 

.523119 

805096 

702717 

E+05 

J 

= 

.220335 

216625 

140614 

E+07 

N  =  15 


U 

= 

.587785 

252292 

473129 

E-00 

U 

= 

.951056 

516295 

153572 

E-00 

V 

= 

.951056 

516295 

153572 

E-00 

V 

= 

.587785 

252292 

473129 

E-00 

A 

= 

-.458173 

664826 

524438 

E  +  10 

A 

= 

. 142370 

102469 

082740 

E+12 

B 

= 

.669809 

170853 

876333 

E+09 

B 

= 

.477908 

042195 

991473 

E+08 

C 

= 

-.378353 

751490 

621970 

E  +  09 

C 

= 

-.193499 

139338 

360761 

E+11 

D 

= 

-.378682 

342664 

795434 

E+08 

D 

= 

-.454582 

176821 

469421 

E+07 

E 

= 

.411125 

234864 

518541 

E+08 

E 

= 

.210710 

209837 

887250 

E+10 

F 

= 

.114632 

583877 

551993 

E  +  07 

F 

= 

.114632 

583877 

551993 

E+07 

H 

= 

-.351218 

040420 

511910 

E+07 

H 

= 

-.164997 

717895 

796311 

E+09 

I 

= 

.125046 

667004 

022639 

E  +  06 

I 

= 

.182440 

628429 

173530 

E+05 

J 

= 

.164584 
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130649 

E  +  06 

J 

= 
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N  -  CUBE  W(X)  =  l/[[[l-(x  )2p 


N  =  16 


u 

= 

•587785 

252292 

473129 

E-00 

U 

= 

.951056 

516295 

153572 

E~00 

V 

= 

.951056 

516295 

153572 

E-00 

V 

= 

.587785 

252292 

473129 

E-00 

A 

= 

-•  169140 

531644 

690936 

E+11 

A 

= 

.605400 

661482 

218181 

E+12 

B 

= 

•244185 

215373 

923036 

E+10 

B 

= 

.187088 

694874 

421396 

E+09 

C 

= 

-.150476 

679570 

194629 

E+10 

C 

= 

-.770393 

426676 

730467 

E+11 

D 

= 

-.128394 

863145 

966847 

E+09 

D 

= 

-.156566 

902027 

263993 

E+08 

E 

= 

.153528 

858854 

910369 

E+09 

E 

= 

.781972 

461899 

526086 

E+10 

F 

= 

.360128 

883371 

733118 

E+07 

F 

= 

.360128 

883371 

733118 

E+07 

H 

= 

-.120846 

820755 

454834 

E  +  08 

H 

= 

-.567722 

820152 

944836 

E+09 

I 

= 

.362626 

791152 

978513 

E  +  06 

I 

= 

.529065 

358144 

058532 

E+05 

J 

.517702 

000258 

462910 

E+06 

J 

= 

.221981 

198739 

585491 

E+08 

N  =  17 


U 

= 

.587785 

252292 

473129 

E-00 

U 

.951056 

516295 

153572 

E-00 

V 

= 

.951056 

516295 

153572 

E-00 

V 

= 

.587785 

252292 

473129 

E-00 

A 

-.613033 

556880 

169686 

E+11 

A 

= 

.251913 

210299 

464150 

E+13 

B 

= 

.881084 

434677 

001421 

E+10 

B 

= 

.715360 

554824 

419001 

E+09 

C 

= 

-.588727 

395465 

418023 

E+10 

C 

= 

-.301372 

432360 

668405 

E+12 

D 

= 

-.432984 

226058 

143313 

E+09 

D 

= 

-.535084 

233106 

935375 

E+08 

E 

= 

.565488 

561952 

925141 

E+09 

E 

= 

.286466 

979652 

642085 

E+11 

F 

s 

.113137 

825434 

613221 

E+08 

F 

= 

.113137 

825434 

613221 

E+08 

H 

= 

-.412664 

656842 

412452 

E+08 

H 

= 

-.193864 

547942 

149012 

E+10 

I 

.105785 

240143 

239198 

E+07 

I 

= 

.154338 

585367 

035380 

E+06 

J 

.162814 

093123 
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